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Abstract 

Let X be the constrained random walk on Zy d G {2,3,4,...}, representing the 
queue lengths of a stable Jackson network and let x G Z^ be its initial position {X is 
a random walk with independent and identically distributed increments except that its 
dynamics are constrained on the boundaries of Z^ so that X remains in Z'}.; stability 
means that X has a nonzero drift pushing it to the origin). Let t„ be the first time 
when the sum of the components of X equals n. The probability = Px{Tn < tq) is 
one of the key performance measures for the queueing system represented by X and its 
analysis/computation received considerable attention over the last several decades. The 
stability of X implies that decays exponentially in n. Currently the only analytic 
method available to approximate Pn is large deviations analysis, which gives the expo¬ 
nential decay rate of Finer approximations are available via rare event simulation. 
The present article develops a new method to approximate Pn and related expectations. 
The method has two steps: 1) with an affine transformation, move the origin to a point 
on the exit boundary associated with r„; let n —> oo to remove some of the constraints 
on the dynamics of the walk; the first step gives a limit unstable / transient constrained 
random walk Y 2) construct a basis of harmonic functions of Y and use this basis to 
apply the classical superposition principle of linear analysis (the basis functions can be 
seen as perturbations of the classical Fourier basis). The basis functions are linear com¬ 
binations of log-linear functions and come from solutions of harmonic systems] these are 
graphs with labeled edges whose vertices represent points on the interior characteristie 
surface H oiY] the edges between the vertices represent conjugacy relations between the 
points on the characteristic surface, the loops (edges from a vertex to itself) represent 
membership in the boundary characteristic surfaces. Characteristic surfaces are alge¬ 
braic varieties determined by the distribution of the unconstrained increments of X and 
the boundaries of Z^}. Each point on T-L defines a harmonic function of the unconsrained 
version of Y. Using our method we derive explicit, simple and almost exact formulas for 
Px^Tn < To) for X representing d-tandem queues, similar to the product form formulas 
for the stationary distribution of X. The same method allows us to approximate the 
Balayage operator mapping / to a: —> Ea, [f{Xx„)l{T^<ro}] Dr a range of stable con¬ 
strained random walks representing the queue lengths of a queueing system with two 
nodes (i.e., d = 2). We provide two convergence theorems; one using the coordinates 
of the limit process and one using the scaled coordinates of the original process. The 
latter is given for two tandem queues (i.e., when the set of possible increments of X is 
{(0,1), (—1,1)(0, —1)}) and uses a sequence of subsolutions of a related Hamilton Jacobi 
Bellman equation on a manifold; the manifold consists of three copies of , the zeroth 
glued to the first along {x : a:(l) = 0} and the first to the second along {x : x{2) = 0}. We 
indicate how the ideas of the paper relate to more general processes and exit boundaries. 
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1 Introduction 


Constrained random walks arise naturally as models of queueing networks and this paper 
treats only walks associated with Jackson networks. But the approach of the paper applies 
more generally, see subsection 19.31 

Let X denote the number of customers in the queues of a d-node Jackson network at ar¬ 
rival and service completion times {Xk{i) is the number of customers waiting in queue i of the 
network right after the arrival/service completion); mathematically, X is a constrained 
random walk on i.e., it has independent increments except that on the boundaries of 

the process is constrained to remain on Z^ (see m for the precise definition of the 
constrained random walk X). Define 

An= IxeZ'l 

I i=i 

and its boundary 

dAn = |x G Z+ : ^ x{i) 

I i=i 

Let Tn be the first time X hits dAn- One of the “exit probabilities” that the title refers to 
is pn = Pxijn < To), the probability that starting from an initial state x G An the number of 
customers in the system reaches n before the system empties. One of our primary aims in 
this paper will be the approximation of this probability. The set An models a systemwide 
shared buffer of size n (for example, if the queueing system models a set of computer pro¬ 
grams running on a computer, the shared buffer may be the computer’s memory) and Tn 
represents the first time this buffer overflows. If we measure time in the number of inde¬ 
pendent cycles that restart each time X hits 0, pn is the probability that the current cycle 
finishes successfully (i.e., without a buffer overflow). 

One can change the domain An to model other buffer structures, e.g., {x g : x{i) ^ n} 
models separate buffers of size n for each queue in the system. The present work focuses on 
the domain An- The basic ideas of the paper apply to other domains, and we comment on 
this in the conclusion. 

For a set a and Ta = {k : X^ G a}, the distribution Ta of X^-^ on a is called the Balayage 
operator. Ta maps bounded measurable functions on a to harmonic functions on T: 

Ta: f ^ g,g{x) = [f{XrJ l{r,<00}] • 

The computation of pn is a special case of the computation of (the image of a given function 
under) the Balayage operator: for a = A^ u dAn u {0}, Ta becomes a tq and if we set 

/ ~ 

{Taf){x) = Ea, [f{X-rJl{r^< 0 D}], X G An, equals Px{Tn < Tq). 

We assume that X is stable, i.e., the total arrival rate t'* is less than the total service 
rate pi for all nodes i of the queueing system that X represents {p and u are linear functions 
of the distribution of the increments of X, see m and m for their definitions). For a 
stable X, the event {r^ < tq} rarely happens and its probability pn decays exponentially 
with buffer size n. The problem of approximating pn has a long history and an extensive 
literature; let us mention two of the main approaches here. The first is large deviations (LD) 
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analysis mmm which gives the exponential decay rate of pn as the value function of a limit 
deterministic optimal control problem (see below). If one would like to obtain more precise 
estimates than what LD analysis gives, the popular method has so far been simulation with 
variance reduction such as importance sampling, see, PQ Chapter VI] and [51[S]; the use of IS 
for similar problems in a single dimension goes back to |16] . The goal of this paper is to offer 
a new alternative, which, in particular, allows to approximate the Balayage operator for a 
wide class of two dimensional systems and gives an almost exact formula for the probability 
Px{Tn < T'o) for tandem networks in any dimension. We explain its elements in the following 
paragraphs. 

One way to think of the LD analysis is as follows, pn itself decays to 0, which is trivial. 
To get a nontrivial limit transform Pn^oVn = \ogpn] using convex duality, one can write 
the — log of an expectation as an optimization problem involving the relative entropy [1] and 
thus Vn can be interpreted as the value function of a discrete time stochastic optimal control 
problem. The LD analysis consists of the law of large numbers limit analysis of this control 
problem; the limit problem is a deterministic optimal control problem whose value function 
satisfies a first order Hamilton Jacobi Bellman equation (see (II48p of Section [7|). Thus, LD 
analysis amounts to the computation of the limit of a convex transformation of the problem. 

We will use another, an affine, transformation of X for the limit analysis. The proposed 
transformation is extremely simple: observe X from the exit boundary. The most natural 
vantage points on the exit boundary dA^ are the corners {ncj, i = 1, 2,3,..., d], where e* are 
the standard basis elements of 

y" = T„(V), r„:M''-M^r„(a;)=y, = (3) 

lx[j) otherwise, 

j = 1,2, - ■■ ,d. Tn is affine and its inverse equals itself. V”, i.e., the process X as observed 
from the corner ne*, is a constrained process on the domain Dy = X (n — Z+) X Z!^; it 
is the same process as X, except that V"' represents the state of the queue not by the 
number of customers waiting in queue i but by the number of spots in the buffer not occupied 
by the customers in queue i. Tn maps the set An to Bn cz Dy, Bn = Tn{An)] the corner ncj 
to the origin of Dy; the exit boundary dAn to dBn = {y ^ y(j)}; finally 

the constraining boundary {z e Zy, z{i) = 0} to 

{y e ZV^ X Z X Z^"* : y{i) = n). 

As n —> 00 the last boundary vanishes and converges to the limit process Y on the domain 
Dy = X Z X Zy“* and the set Bn to 


B = 


|y e ^ 


t yU) 


( 4 ) 


Figured] sketches these transformations for the case of X representing lengths of two tandem 
queues and for i = 1 (the random walk X represents tandem queues if its set of possible 
jumps are ei, —e* + Cj+i , i = I, 2,3,..., d — 1 and —e^, if X is of this form we will call it a 
“tandem walk;” for the exact definition, see (fT^ L 
The boundary of B is 

= jy 6 Dy,y(i) = ^ (5) 

V j = J 
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Figure 1: The transformation 


dB 
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the limit stopping time 

r = inf{A: : Yfc e dB] (6) 

is the first time Y hits dB. The stability of X and the vanishing of the boundary constraint 
on i implies that Y is unstable / transient, i.e., with probability 1 it wanders off to oo. 
Therefore, in our formulation, the limit process is an unstable constrained random walk in 
the same space and time scale as the original process but with less number of constraints. 

Fix an initial point y e B m the new coordinates; our first convergence result is Propo¬ 
sition [3T] which says 

Pn = Pxffrn < To) ^ By {t < CO), (7) 

where Xn = Tffy). The proof uses the law of large numbers and LD lowerbounds to show 
that the difference between the two sides of © vanishes with n. With © we see that the 
limit problem in our formulation is to compute the hitting probability of the unstable Y to 
the boundary dB. 

The convergence statement ([7|) involves a fixed initial condition for the process Y. In 
classical LD analysis, one specifies the initial point in scaled coordinates as follows: Xn = 
]nx\ e An for x e Then the initial condition for the Y^ process will be = Tffxn) (thus 
we fix not the y coordinate but the scaled x coordinate). When Xn is defined in this way, 
(j7|) becomes a trivial statement because its both sides decay to 0. For this reason, Section [3 
studies the relative error 

\PxniTn < To) - Pyffr < CO)\ 

Pxn{rn<ro) ’ 

Proposition 17.11 says that this error converges exponentially to 0 for the case of two dimen¬ 
sional tandem walk (i.e., the process X shown in Figured]). The proof rests on showing 
that the probability of the intersection of the events {r^ < tq} and {r < oo} dominate the 
probabilities of both as n ^ oo. For this we calculate bounds in Proposition 17.31 on the LD 
decay rates of the probability of the differences between these events using a sequence of 
subsolutions of a Hamilton Jacobi Bellman equation on a manifold', the manifold consists of 
three copies of 1R+, zeroth copy glued to the first along di, and the first to the second along 
^ 2 , where di = {x e : x{i) = 0}. Extension of this argument to more complex processes 
and domains remains for future work. 

For a process X in d dimensions, each affine transformation Tf, i = 1, 2, 3,..., d, gives a 
possible approximation of Px{Tn < To)- A key question is: which of these best approximates 
Pxijn < To) for a given xl Proposition 17.11 savs that, for the two dimensional tandem walk, 
i = 1 works well for all points x = ]nx\ as long as x(l) >0. In general this will not 
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be the case (i.e., depending on x, removing one constraint of the process may give better 
approximations than removing another); subsection 17.11 comments on this problem. 

The convergence results ([7|) and ([8]) reduce the problem of calculation of Px{Tn < To) to 
that of Pyij < oo). This constitutes the first step of our analysis and we expect it to apply 
more generally; see subsection 19.31 Computation of Py{T < oo) is a static linear problem and 
can be attacked with a range of ideas and methods. 

Sections H El and [6] apply the principle of superposition of classical linear analysis to 
the computation of Py{T < oo) and related expectations. The key for its application is to 
construct the right class of efficiently computable basis functions to be superposed. The 
construction of our basis functions goes as follows: the distribution of the increments of 
Y is used to define the charaeteristic polynomial p : C” —> C. p can be represented both 
as a rational function and as a polynomial; to simplify our analysis we use the polynomial 
representation in two dimensions and the rational one in d dimensions. In the rest of this 
paragraph we will only refer to the higher dimensional definitions; the definitions for the case 
of two dimensions are given in subsection 14.11 We call the 1 level set of p, the charaeteristic 
surface of Y and denote it with T-L, see (j98p . T-L is, more precisely, ad — 1 dimensional complex 
afhne algebraic variety of degree d + 1. Each point on the characteristic surface P defines 
a log-linear function (see (195 h i that satisfies the interior harmonicity condition of Y (i.e., 
defines a harmonic function of the completely unconstrained version of T); similarly, each 
boundary of the state space of Y has an associated characteristic polynomial and surface, p 
can be written as a second order polynomial in each of its arguments; this implies that most 
points on P come in conjugate pairs, there are d — 1 different conjugacy relations, one for 
each constraining boundary of Y. The keystone of the approach developed in these sections 
is the following observation: log-linear funetions defined by two points on P satisfying a 
given type of eonjugacy relation ean he linearly combined to get nontrivial funetions whieh 
satisfy the corresponding boundary harmonieity eondition (as well as the interior one); see 
Proposition 15.11 Based on this observation we introduce the concept of a harmonie system 
(Definition 15.2p which is an edge-eomplete graph with labeled edges representing a system 
of variables and equations: the vertices represent the variables constrained to be on P, the 
edges between distinct vertices represent the conjugacy relations between the variables that 
the edges connect (the label of the edge determines the type of the conjugacy relation) and 
its loops (an edge from a vertex to itself) represent membership on a boundary characteristic 
surface (the label of the loop determines which boundary characteristic surface). We show 
that any solution to a harmonic system gives a harmonic function for Y in the form of 
linear combinations of log-linear functions (each vertex defines a log-linear function). The 
computational complexity of the evaluation of the resulting harmonic function is essentially 
determined by the size of the graph. 

In two dimensions (Section [3|) edge-complete graphs have 1 or 2 vertices and the above 
construction gives a rich enough basis of harmonic functions of Y to approximate the image 
of any function on dB under the Balayage operator; with the use of these basis functions 
the approximation of Ej^[/(W)l{.r<oo}] for any given bounded / reduces to the solution of 
a linear equation in K dimensions, where K is the number of basis functions used in the 
approximation (Section 18.21 gives an example with K = 12). Once the approximation is 
computed, the error made in the approximation is simple to bound when / is constant 
outside of a bounded support and satisfies / > 0 or / < 0. The restrictions of the basis 
functions on dB are perturbed versions of the restriction of the ordinary Fourier basis on Z 
to Z+; for this reason we call the constructed basis a “perturbed” Fourier basis. 

Section [S] gives the definition of a harmonic system for d-dimensional constrained walks 
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and prove that any solution to a harmonic system defines a y-harmonic function (Proposition 
15.2|) . In Section [B] we compute explicit solutions to a particular class of harmonic systems for 
the d dimensional tandem walk; the span of the class contains Py{T < oo) exactly. Hence we 
obtain explicit formulas, similar to the product form formulas for the stationary distribution 
of Jackson networks [TO], for Py{T < oo) in the case of constrained walks representing tandem 
queues in arbitrary dimension (see Proposition 16.51 the two dimensional version of the same 
formula is written out more explicitly in display (|163p i. If we take exponentiation and 
algebraic operations to be atomic, the complexity of evaluating the formula is independent 
of y (and hence of the buffer size n if Pyij < oo) is used to approximate Px{jn < tq) ) and 
depends only on the dimension d ol X. 

Section [8] gives example computations using the approach developed in the paper. Three 
examples are considered: the tandem walk in two dimensions, a non tandem walk in two 
dimensions, tandem walk in 4 and 14 dimensions. The conclusion (Section[B|) discusses several 
directions for future research. Among these are the application of the approach of the present 
paper to constrained diffusion processes (subsection 19.11) , the study of nonlinear perturbed 
second order HJB equations which arise when one would like to sharpen large deviations 
estimates (subsection 19.2p and the sizes of boundary layers which arise in subsolution based 
IS algorithms applied to constrained random walks (subsection 19.51) . 

2 Definitions 

This section sets the notation of the paper, defines the domains, the processes and the 
stopping times we will study and states some elementary facts about them. 

We will denote components of a vector using parentheses, e.g., for x e Z‘1, x(i), i = 
1, 2,3,..., d denotes the component of x. 

For two sets a and 6, denotes the set of functions from b to a. For a function / on 
a set a, we will write f\c to mean f’s restriction to a subset c a. For a finite set a, |a| 
denotes the number of its elements. We will assume that elements of sets are written in a 
certain order and we will index sets as we index vectors, e.g., a(l) denotes the first element 
of a and a(|a|) the last. 

Our analysis will involve several types of boundaries: the coordinate hyperplanes of 
the constraining boundaries of constrained processes and the boundaries of exit sets. To 
keep our notation short and manageable, we will make use of the symbol d to indicate that 
a set is a boundary of some type. 

Define A/q = {0,1, 2,..., d} and J\f+ = A/q ~ {0}; A/V is the set of nodes of X. For a cz A/V 
the coordinate hyperplanes of are 

da = e Z'^, z{j) = 0 Vj G a j . 

We will use the letter a to denote hitting times to these sets; for any process P on Z'^ define 

= inf{A: : Pk e da}; (9) 

in what follows the process P will always be clear from context and we will omit the 
superscript of a. If a = {j} for some j e A/+, we will write aj rather than cryj; the same 
convention applies to dyj. 

We will denote the domain of a process P by H-p; dH-p will denote its constraining 
boundary if it has one; Dp will denote Dp — dHp. 
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We will often express constraints of constrained random walks using the constraining 
maps TTa, a cz A/V, defined as follows: 


'^a{x,v) 


X + V, if x{j) + v{j) ^0 Vj e 
X, otherwise, 


where x 6 and x 6 Z'^. If a = 0, we will write tt instead of tt 0 ', tt constrains to Z^^ any 
process to which it is applied (see the dehnition of X and Y below to see how TTa is used to 
dehne constrained random walks). Other than a = 0 this paper will only use vrpj, i e M+ 
and for most of the paper we will assume i = 1. To ease notation we will write vrj rather 
than TTpj; vTj constrains a given process to be positive in its M+ — {i} coordinates. 

If P is a random walk with increments V{V), constrained to stay in Q-p cz Z'^, and 
S cz 0-p, dehne 

5° = {s e 5 : Up n (s + V{V)) ^ S},dS = S - (10) 


The notation doesn’t state explicitly the T’-dependence of these terms; but whenever we 
use them below, the underlying process V will always be clear from context. In what follows 
V will be either X, T”, Y or Z, all of which are dehned below. 

We want to compute certain probabilities/ expectations associated with a constrained 
random walk. This will involve three transformations of the original process: an affine change 
of variables, taking a limit (this will drop one of the boundary constraints of the process) and 
removing all constraints which makes the process an ordinary random walk with independent 
and identically distributed (iid) increments. We will show the original process with X, the 
result of the affine transformation with T”, the limit process with Y and the completely 
unconstrained process with Z. 

X will denote a constrained random walk on fix = Z'/ with independent increments 
Ik ^ {ci — ej,i Y j ^ A/q} where eo is the zero vector in Z'’* and e* e Z'’*, i ^ 0 is the unit 
vector in the direction i. To keep X in its domain, the increments are constrained on the 
boundaries of Z!^: 


TT 


Xo = X e Z+, 

Xk+l = Xk + 7T{Xk,Ik) 

Ci-ej, if x(j) > 0 
0, otherwise. 


(x,ei - ej) = 


( 11 ) 


where, by convention, “x(0)” means “1” (or some other positive quantity). The constraining 
boundary of X is dUx = {^jeAf+^j) ■ 

We denote the common distribution of the increments Ik by the matrix p, i.e., p{i,j) is 
the probability that Ik equals ei — ej, i Y j ^ A/q; p{i,i) = 0 for i e A/q. V(X) will denote 
the set of increments of Ik with nonzero probability: 


V{X) = {ei - ej : p{i0) > 0}. 

Remark 1. For v = Ci — ej, i,j e Aq x A/q, ^‘p{vy’ will denote p{i,j). 

For our probability space we will take V(X)^. {Ik} are the coordinate maps on V(X)^, 
^ is the (T-algebra generated by {Ik} and P is the product measure on V(X)^ under which 
Ik are an iid sequence with common distribution p. 
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X describes the dynamics of the number of customers in the queues of a Jackson net¬ 
work, i.e., a queueing system with exponentially and identically distributed and independent 
interarrival and service times. In this interpretation, 0 e A/q represents the outside of the 
queueing system and for i e A/+, Xk{i) is the number of customers in queue i at the 
jump of the system (an arrival or a service completion). X is on the boundary di when the 
queue is empty and the constrained dynamics on the boundary means that the server at 
node i cannot serve when its queue is empty. The increment Ci — ej, i,j ^ 0, i ^ j, represents 
a customer leaving queue i after service completion at node i and joining queue j] Ci, i ¥= 0, 
represents an arrival from outside to queue i and —ej a customer leaving the system after a 
server completion at queue j. We will assume that the Markov chain defined by the matrix 
p on A/q is irreducible. One can also represent p as arrival, service and routing probabilities: 

Aj=p(0,j), pj = p{j,k),j e J\f+, r{i,j) = p{i,j)/pi. (12) 

jeXo 


The irreducibility of p implies that 


u = X + ru (13) 

has a unique solution, is the total arrival rate to node i when system is in equilibrium. 
We assume that the network corresponding to X is stable, i.e, 

Pi = -<1, ieN+. (14) 

Pi 


We will pay particular attention to tandem networks, i.e., a number of queues in tandem; 
these are Jackson networks whose p matrix is of the form 

p(0,1) > 0,p(0,j) = OJ ^ l,p(d,0) = pd,p{j,j + 1) = PjJ eM) - {0,d}; (15) 

for tandem queues A will denote the only nonzero arrival rate p(0,1); then = p(0,1) = A 
for all i e A/+. We will call the random walk X a tandem walk if it represents the queue 
lengths of a tandem network. 

The domain is defined as in ([T]). We will assume 

Xj > 0 at least for some j e J\f+. (16) 

With this and (fTOl) . dAn indeed equals the right side of (l2|). Define the stopping times 

Tn = {A: . X]^ G dAji} 


and 

Pn = P{Tn < To). 

If (fT6l) fails, Pn becomes trivial. 

Define the input/output ratio of the system as 

p(a 0) ■ 


( 17 ) 


Stability of X implies 
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Proposition 2.1. 


r < 1. 


(18) 


Proof. By definition 

Yi r{k,j)i'k, 

for j e Af+. Sum both sides over j: 

2 M0,j) = Y Y ^(pk)) 

jeN+ jej\f+ ke{jY 

- S - S 

jeAf+ jeAf+ 


Then 

_ IjjeA4 p{j, 0) 

pip 0 ) ” pip 0 ) 

which is an average of the utilization rates pi which are by assumption all less than 1; (1181) 
follows. □ 

Tn and are defined as in Q. depends on i; when we need to make this dependence 
explicit we will write T^; for most of the analysis of the paper i will be fixed and therefore 
can be assumed constant, and, unless otherwise noted, in the following sections we will take 
i = 1. 

Define 


li 6 Ii{j,k) =0,j ^ k, Ii{j,j) = l,i ^ j, = -1. (19) 

li is the identity operator except that its diagonal term is —1 rather than 1. Then 

Tn = nei+Ii. ( 20 ) 

Define the sequence of transformed increments 

Jk=^iih)\ (21) 

Jfc and Ik take the same values except that Jk = Qj + &i whenever Ik = ej — and Jk = — 
whenever Ik = Ci. Define 


V{Y) = {Z,v,veV{X)}. 

Remark 2. For v e V(y) we will shorten p{Tiv) to p{v) (remember, per Remark [1] for 
ei - ej, i,j G Mo, p{v) denotes p{i,j)). 

The limit unstable constrained process Y is 

Yq = y e Z^, Yfc+i = Yk + T^ifYk, -Ik). (22) 

Y has the same dynamics as Y'^ except that Y has no constraining boundary on its 
coordinate; therefore its state space is Dy = x Z x Z'^~*. The domain B for Y is defined 
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as in (HI) and its boundary dB is defined from B using (IlOp and coincides with the right side 
of (O. Let T be as in ([S]). Define 

Cn = inf I: yfc(f) = 2 YfcO') + n| ; (23) 


note T = Co- 

X, Y"' and Y are all defined on the same probability space (V(X)^,^, P); the measure 
P and their initial positions determine their distributions. We will use a subscript on P to 
denote the initial positions, e.g., Px{jn < tq) is the same as P{Tn < tq) with Xq = x and 
Pyir < oo) means P{t < oo) with Yq = y. 

We note a basic fact about Y here: 

Proposition 2.2. For y 6 « 

PyiCn A Co = «)) = 0. (24) 


Proof. Set 


For y s Z'l, y : y{i) < n, 


c = 2 

jeXo 


PyiCn A Co ^ n) > c = c” > 0 (25) 

because, at least the sample paths whose increments consist only of {—Cj, i e ACf ,p(0, i) > 0} 
push Y to dB in n steps and the probability of this event is c"'. 

hfc = PnfcACnACo is Markov on = jy e ^ Th=i viC) ^ nj (because Y is Markov). 

The boundary of the last set is 

f d d 

SBn = t 2 /e Z+,0 = ^ yii) or ^ yii) = n 

( i=l i=l 

By dehnition 

PyiCo Cn = ^ PyiYk^ Bn- dBn)- (26) 

The bound (f2^ implies Py e Bn — dBr^ ^ 1 — c'. This and that Y is Markov give 
Py e Bn — dBn^ ^ (1 — c')^. This and (f26]l imply 

^’j/(CoACn = a))^(l-c')". (27) 

Letting A: —> oo gives ((Ml) . □ 

3 Convergence - initial condition set for Y 

This section shows that the affine transformation of observing the process from the exit 
boundary really gives approximations of the exit probabilities we seek to compute. The 
present convergence result specifies the initial point for the Y process. This allows a simple 
argument that works for general stable X and uses LD results only roughly to prove that 
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certain probabilities decay to 0. In Section [7] we will prove a second convergence result (for 
the case of two dimensional tandem walk) where the initial point is given for the X process; 
this will require a finer use of large deviations decay rates. 

Denote by X the law of large numbers limit of X , i.e., the deterministic function which 
satisfies 

limP^„ ( sup \Xk/n - > ,5 ) = 0 (28) 

^ \k^ton J 

for any <5 > 0 and to > 0 where Xn £ is a sequence of initial positions satisfying ^ —> 
X G (see, e.g., [Ml Proposition 9.5] or [21 Theorem 7.23]). The limit process starts from 
Xq = X, is piecewise affine and takes values in then st = Yii=i starts from x(f) is 
also piecewise linear and continuous (and therefore differentiable except for a finite number 
of points) with values in M_|_. The stability and bounded iid increments of X imply that s is 
strictly decreasing and 

Cl > —s > Co > 0 (29) 

for two constants ci and cq. These imply that X goes in finite time to 0 g and remains 
there afterward. 

Fix an initial point y G Dy for the process Y and set Xn = Tn{y)] dTH]) implies 



n 


(30) 


Proposition 3.1. Let y and be as above. Then 


lim Px„('rn < To) = Py{T < oo). 

n —>00 

Proof. Note that for n > y{i), Xn g A^- Define 

Mfc = maxl)(i), = minX;(z). 

l^k l^k 

M is an increasing process and Mr is the greatest that the component of Y gets before 
hitting dB (if this happens in finite time). The monotone convergence theorem implies 

P„(t < oo) = lim Py{T < 00, Mr < n). 
nyoo 


Thus 


Py{T < oo) = Py{T < 00, Mr < n) + Py{T < 00, Mr ^ u) 
and the second term goes to 0 with n. Decompose Px^i^n < tq) similarly using M^: 

PxS^n < To) = Px^ {rn < Tq, > O) + Px^ (t^ < Tq, = O) . 


(31) 


On the set {M^ > 0}, the process X cannot reach the boundary di before r^, therefore over 
this set 1) the events {t„ < tq} and {r < oo} coincide (remember that X and Y are defined 
on the same probability space) 2) the distribution of {Tn{X),n — M^) is the same as that 
of {Y,M) upto time r^. Therefore, 

= Py{T < 00, XIr <n) + Pxr (Tn < Tq, = O) . 
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The first term on the right equals the first term on the right side of (|3ip . We know that the 
second term in (I3ip goes to 0 with n. Then to finish our proof, it suffices to show 

lim (r„ < To, = o) = 0. (32) 

0 means that X has hit di before r^. Then the last probability equals 

Px^ (o-i <Tn< To) , (33) 

which, we will now argue, goes to 0 (ui is the first time X hits di\ see ([9|)). (i30p implies 
To = ej. Define T = inf{f : Xt{i) = 0} and = inf{f ; = 0 e By definition T ^ < oo 

Now choose to in (l28p to be equal to t®, define = {sup/j^jo„ e \Xk/n — > <5} and 

partition (l33P with 

Pxr, (c^i <Tn< To) = Pxn <Tn< Tq} Pi Cn) + Pxr, < Tn < To} P C^) ■ (34) 

The first of these goes to 0 by (I28p . The event in the second term is the following: X remains 
at most nS distance away nX until its nt^ step, hits di then dAn and then 0. These and (I29p 
imply that, for n large enough, any sample path lying in this event can hit dA^ only after 
time nt^. Thus, the second probability on the right side of (|34p is bounded above by 

PxA{nt^ <Tn< To} pC^). 

The Markov property of X, {cjj < Tn < tq} cz {r„ < To} and (l28p imply that the last 
probability is less than 

Px{Tn < '^o)Pxn{^nfi ~ ^)- 

x:\x\^n5 

For |x| ^ n6, the probability Pxijn < To) decays exponentially in n [HI Theorem 2.3]; then, 
the above sum goes to 0. This establishes (13211 and finishes the proof of the proposition. 

□ 


4 Analysis of y, d = 2 

Let us begin with several definitions for the general dimension d; because we will almost 
exclusively work with the Y process from here on, we will shorten V(y) to V. A function 
V ; ^ C is said to be a harmonic function of the process Y (or Y-harmonic) on a set 

O c Dy if 

y{y) = [^(y)] = Yi ^ o, (35) 

veV 

where we use the convention set in Remark [2j Throughout the paper O will be either B° 
or Dy, and the choice will always be clear from context; for this reason we will often write 
“...is T-harmonic” without specifying the set O. V is said to be Z-harmonic on O c Z*^ if 

V(z) = E, [l/(Zi)] = 2 ^(2 + v))p(v),z e O. 

veV 

Z-harmonicity and R-harmonicity coincide on Qy. 

Above we have assumed the domain of V to be Z*^. If V is defined only on a subset 
a c Z'^, it can be trivially extended to all Z*^ by setting it to 0 on Z'^ — a. Thus, the above 
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definitions can be applied to any function defined on any subset of we will use a similar 
convention for most of the definitions below. 

The Markov property of Y implies that 

h:y^Ey[f(Yr)l(r<oo}],y^S, (36) 

is a harmonic function of Y whenever the right side is well defined for all y e B°. Note that 
h is the image of the function / under the Balayage operator The dynamics of Y 

and the definition of B imply that is a distribution on dB\ for this reason we will call 

harmonic functions of the form ([36]) 5i?-determined. If a function / is defined over a domain 
larger than B, we will call / 5i?-determined, if its restriction to B is so. 

The analysis of the previous section suggests that we approximate 

PxiYn ^ T)) 

with W{Tn{x)) where 

W(y) = Py{T <co) (37) 

for any stable Jackson network X. IT is a dB determined harmonic function of Y, in 
particular it solves ([35]l with O = B°. That W{y) = Py{T < oo) = 1 for y e dB implies that 
W also satisfies the boundary condition 


V\dB = 1- (38) 

Then IT is a solution of ()35l38p with O = B° . Large deviations analysis of IT is an 
asymptotic analysis of the system (I35I38P that scales V to log V and uses a law of large 
numbers scaling for space and time. With the y coordinates, we no longer need to scale V, 
time or space and can directly attempt to solve (I35I38I) - perhaps approximately. We have 
assumed that X is stable; this implies that T-Afc; k = 1,2, 3,..., is unstable and therefore, the 
Martin boundary of this process has points at infinity. Then one cannot expect all harmonic 
functions of Y to be 5i?-determined and in particular the system (|35l38p will not have a 
unique solution; hence, once we find a solution of (l35l 1381) that we believe (approximately) 
equal to Py{T < oo), we will have to prove that it is 5i?-determined. 

Define 

Bz = \zeZ'^: z{l) ^ j] z{j) 

I 1=2 

The unconstrained version of (135 p is 

V{z)=E,[V{Zi)] = Y,V{z + v)p{v), (39) 

veV 

z e O and that of ([5SP is 

V\8B, = 1. (40) 

A function is said to be a harmonic function of the unconstrained random walk Z on O if it 
satisfies (l39p . 

Introduce also the boundary condition 

V\sB = f, f:dB^C, (41) 

for Y which generalizes ()38h . 

Our idea to [approximately] solve (|35l38p is this: 
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1. Construct a class Ty of “simple” harmonic functions for the process Y (i.e., a class of 
solutions to For this 

(a) Construct a class Tz of harmonic functions for the unconstrained process Z (i.e., 
a class of solutions to (j39p with O = Z"^), 

(b) Use linear combinations of elements of Tz to find solutions to (13511 . 

2. Represent [or approximate] the boundary condition (|38p by linear combinations of the 
boundary values of the dR-determined members of the class Ty. 

The definition of the class Tz is given in (I5TI1 and that of Ty is given in (IHUIl . 

This section treats the case of two dimensions, where this program yields, for a wide class 
of processes, approximate solutions of (l3^ not just with / = 1, i.e, the boundary condition 
(i38l) . but with any bounded /, i.e., the boundary condition (HTl) . 

Stability of X implies p(2, 0) a p(l, 0) > 0 and we will assume 

p(2,0)>0; (42) 

otherwise one can switch the labels of the nodes to call 2 the node for which p(z, 0) > 0. 

4.1 The characteristic polynomial and surface 

Let us call 

p(/3,a) = e C^, (43) 

veV 

the characteristic polynomial of the process Z for Bz, 

p(/3, a) -/3a = 0 (44) 

the characteristic equation of Z for Bz and 

% = {(/3,a) : p(/3,a) — /3a = 0} 

the characteristic surface of Z for Bz- We borrow the adjective “characteristic” from the 
classical theory of linear ordinary differential equations; the development below parallels that 
theory. Figure [2] depicts the real section of the characteristic surface of the walk whose p 
matrix equals 

/ 0 0.05 0.1 \ 

p = 0.35 0 0.12 . (45) 

\0.3 0.08 0 / 

B is an affine algebraic curve of degree 3 [U Definition 8.1, page 32]; its d dimensional 
version in Section [5] will be an affine algebraic variety of degree d + 1. We will need, for the 
purposes of the present paper, only that points on these varieties come in conjugate pairs 
(see below). A thorough study/ description of the geometry of these varieties (and their 
projective counterparts) and its implications for constrained random walks will have to be 
taken up in future work. 

p is a second order polynomial in a [/3] with second and hrst order coefficients in /3 [a]: 
p(/3, a) = (p(l, 0)a + p{2, 0))/3^ + (p(l, 2)a^ + p{2, 1) - a)/3 + (p(0, 2)a^ + p(0, l)a) (46) 

p(/3, a) = (p(0,2) + /3p(l, 2))a^ + (p(0,1) + p(l, 0)/3^ - /3)a + {p{2, 0)/3^ + p{2, l)/3). (47) 
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A singularity For /3p(l,2) +p{0,2) = 0 (HTt) becomes affine. Ifp(l,2) = p{0,2) = 0, (H71) 
is affine for all values of (3 and the method developed below is not applicable. But such walks 
are essentially one dimensional (only their first component can freely move and their second 
component decreases to 0 and stay there upon hitting it) and yield to simpler methods. In 
what follows the /3 we will work with will always satisfy 

/3p(l,2)+p(0,2) ^0. 



Figure 2: The real section of the characteristic surface T-i of the walk defined by p of (14511 : 
the end points of the dashed line are two conjugate points, see dZZ]) 


4.2 log-linear harmonic fnnctions of Z 


The Z-version of the random times Cn of (123 p and r of ([6]) are 

= inf {k : Zk 5 dBz} 

Cn = inf I A: : Zk{i) = ^ Zk{j) + n 

I 



We will omit the Z superscript below because the underlying process will always be clear 
from context. 

For T < CO, Zr takes values in dBz and Zt(1) = Zt-(2). Therefore, the distribution of Z^ 
on dBz is equivalent to the distribution of Zt(1) on Z whose characteristic function is 




{t<oo} 

That Zt(1) is integer valued makes the above characteristic function periodic with period 
27r therefore we can restrict 0 e [0, 27r); setting a = we rewrite the last display as 


a 




{r<oo} 


, Q; S 


s\ 


(48) 


where = {ueC: |u| = 1} is the unit circle in C. For each fixed a e the right side of 
()48ll dehnes a harmonic function of the process Z on B^ as 2 : varies in this set. Our collection 
of harmonic functions Bz for the process Z will consist of these and its generalizations when 
we allow a to vary in C. For a e C the function 2 ; —> is an eigenfunction of the 

translation operator on and Z is a random walk on the same group. These imply 
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Proposition 4.1. Suppose 
for a e C. Then 
for z e C where 


E, 


( 1 , 0 ) 


< 00 


*^1 ■‘-{t<oo} 


( 49 ) 


(7 = E, 


( 1 , 0 ) 


1{t<00} 


Furthermore, {U,a) is on the characteristic surface TL. 

Proof. The proof will be by induction on z{l)—z{2). ()4U]1 is true by definition for z(l)—z{2) = 
0. Assume now that holds for z{l) — z{2) = fc — 1 ^ 0 and fix z with z{l) — z{2) = k. 
The invariance of Z under translations implies 


= Z + {j - 1, j),Cfc-l < oo) = P(i,o)(^r = {jJ),T < oo) 


(50) 


for j e Z. The strong Markov property of Z and Cfc-i < i" imply 



= E^ 

l{Cfc-l<C0}E2 

Cl ^l{r<co}^^4_i 


= E^ 

1{4_i<co}E2 

Cl ^l{r<co} 


The random variable is discrete; then, one can write the last expectation explicitly as 

the sum 


Xj ®'2+0'-i,j) 

j=-co 


Is 


= z + {j . 


z' = z {j — 1, j) satisfies ^'(1) — z'{2) = k — 1] this, the induction hypothesis and ([SCT]) give 


00 

j^-CO 

00 

= (j,j),l{r<oo}) 

j=-co 

By definition, the last sum equals E(i o) l{r<oo}] = U and therefore 

= C/^(^)“^“^(2)Q,^(2)fy ^ jjz{l)-z{2) ^z{2) ^ 

i.e., (|49p holds also for 2 with z{l) — z{2) = k. This finishes the induction and the proof of 
the first part of the proposition. 

The Markov property of Z and the first part of the proposition imply that g : z ^ 
[/ 2 (i)- 2 ( 2 )q, 2 ( 2 ) -g ^ harmonic function of Z on i.e., it satisfies (I39p . Substituting g in 
([39]) implies that {U,a) is on the characteristic surface TL. □ 

Conversely, any point on TL defines a harmonic function of Z: 

Proposition 4.2. For any (/3,a) e TL, z ^ ^ ^ ^ harmonic function 

of Z. 
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Proof. Condition Z on its first step and use p(/3, a) = 13a. 
For {(3, a) e C^, define 

[(/3,a),-] : ^ C, [(/3,a),z] = 

The last proposition gives us the class of harmonic functions 

Pz = {[W,a),-], {f3,a)en} 


□ 


(51) 


for Z. 

4.3 ^S^-determined harmonic functions of Z 

A harmonic function h of Z on Bf. is said to be di?^-determined if 

h{z) = Ex [/ {Zr) l{r<oo}] ) 2: e 

for some / : dBz —> C for which the right side is well defined for all z e Bf. The above 
display defines the Balayage operator of ^ on dBz, mapping / to h; thus, h is dBz- 

determined if and only if it is the image of a function / under the Balayage operator 7(5° )c. 
In our analysis of Y and its harmonic functions we will find it useful to be able to differentiate 
between harmonic functions of Z which are dB^-determined, and those which are not. The 
reader can skip this subsection for now and can return to it when we refer to its results in 
subsection liTl 

For each a e C satisfying 

ap{l,0)+p{2,0) AO (52) 

(m is a second order polynomial equation in f3 (see (1461) 1 with roots 

^ a - p{f, 2)a^ - p(2, 1) - a/A ^ a - p(l, 2)a^ - p(2, 1) + 

2(p(2,0)+p(l,0)a) ’ 2(p(2,0)+p(l,0)a) ’ ^ ^ 

where 

A(a) = (p(l, 2)a^ + p(2,1) - if - 4(p(2, 0) + p(l, 0)a)(p(0, l)a + p(0, 2)a^), 
and ^/z denotes the complex number with nonnegative real part whose square equals z. 

Remark 3. Unless otherwise noted, we will assume (1521) . The stability condition (|14p rules 
out p(l, 0) = p{2, 0) = 0; if p{l, 0) = 0, (f^ always holds. If p{l, 0) ^ 0 and p{2, 0) = 0, (1^ 
fails exactly when a equals 0, a value which represents a trivial situation (Balayage of the 
zero function on dBz). When p(l,0),p(2,0) A 0, ([5^ fails only for a = < 0. This 

value may be of interest to us in the next subsection and we comment on it there in Remark 

El 


Our next step is to identify a set of a’s for which one of the roots above gives a dBz- 
determined harmonic function. In this, we will use m Exercise 2.12, Chapter 2, page 54] 
(rewritten for the present setup): 

Proposition 4.3. For a function f : dBz —> M+ z —>■ [f{Zr)l{r<oo}] > ^ ^ Bz, is the 
smallest function equal to f on dBz and harmonic on B'^. 
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We begin with a = 1. 


Proposition 4.4. 






{t<oo} 


= P,(r<(») =/3i(l)^W =r^«, 


where r is the input/output ratio (fT71) of X . 

Proof. 

A(l) = {{p{l, 0) + p{2, 0)) - (p(0,1) + p{0, 2)))2 
Proposition 12.11 implies 

= (p(l, 0 ) + p{2, 0 )) - (p( 0 , 1 ) + p( 0 , 2 )) > 0 . 


(54) 


(55) 


Then /3i = r < f32 = 1 and (l,r) and (1,1) are points on the characteristic curve Ti. 
Proposition 14.11 implies P(i^o)('^ < oo) = /3i = r or P(i^o)('^ < oo) = /32 = 1- Proposition 14.31 
implies that the former must hold. Proposition 14.11 now implies (|54p . □ 

For a complex number z let 3ft(z) [9(z)] denote its real (imaginary) part. If we write 
a = e (0, 27r) and set x = cos{9) then 

5P(A) = 2x2(p(2, ^^2 ^ 2)2) - 2x(p(2,1) + p(l, 2) + 2p(l, 0)p(0, 2) + 2p{2, 0)p(0,1)) 

+ 1 - 4(p(0, l)p(l, 0) + p{2,0)p{0, 2)) - (p(2,1) - p{l, 2))2 


9(A) = sin(0)(2x(p(l, 2)2 - p{2, 1)2) + 4p(0, l)p(2,0) - 4p(0, 2)p(l, 0) + 2p(2,1) - 2p(l, 2)). 

9(A)/sin(0) is affine in x; to simplify exposition, we will assume that this function has a 
unique root lying outside of the interval (— 1 , 1 ): 

2j.(0,2)p(l,0)-2j.(0,l)p(2,0)+p(l,2)-p(2,l) ,,, 

p(l, 2 P-p( 2 ,l )2 <“> 

See the end of this subsection for comments on ()56H . This and sin(0) 7 ^ 0 imply 

9(A(a)) 7 ^ 0,6» 6 ( 0 , 7 r). (57) 


Proposition 4.5. 


/ 3 i{a) ¥= /32(a) 


for a = e*®. 

Proof fi- 132 = - p( 2 , 0 )^( 1 ,o)a and ([58]) will follow from 


A(e*®) 7 ^ 0,0 e [0,27r). 


(58) 


(59) 


A is a polynomial with real coefficients. Then A(exp(—/0)) = A(exp(z0)) and hence, it 
suffices to prove (l5^ for 6 e [ 0 , 7 r]. ([57)1 implies ([5^ for 6 e [ 0 , 7 r). For 0 = 0: a = 1, (1^ 
and Proposition 12.11 imply A(l) > 0. For 0 = tt: a = —1 and 


A(-1)>A(1). 


These prove (l5^ for 0 6 [0, vr] and complete the proof. 


( 60 ) 

□ 
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Proposition 4.6. 


(61) 


9^(5i (e*®) , 0 6 (0, 27r] 9 ^ ^2 (e*®) , 9 e {Q, 27r] 


are continuous. 

Proof. We have defined to mean the complex number with positive real part whose square 
equals z; this definition leads to a discontinuity only when z passes the negative side of the 
real axis on the complex plane. Then the only possibility of discontinuity for the functions 
/3i and /32 (as functions of 9, as in (I6ip l is if A(exp(f0)) crosses this half line as 9 varies in 
[0, 27r). But (1551) . (l57|) and (l6n|) imply that as 9 varies in [0, vr), A defines a curve C starting 
from and ending at the positive real line and lying on either on the positive or the negative 
complex half plane. That A is a polynomial with real coefficients implies that A(exp(i0)), 
9 e [vr, 27r) is the mirror image C of C with respect to the real line; C and C together dehne a 
closed loop that crosses the real line twice on its positive side. These imply that VA defines 
a continuous closed loop in {z e C : 9?(z) > 0}, from which the statement of the theorem 
follows. □ 

Figure [3] depicts the curves traced on the C-plane by /3i, -s/A and /I 2 as 9 varies in [0, 27r) 
for the p matrix of ()45h . 



Figure 3: /3i(e®®), ^/and a section of /32(e*^) (intersecting ^/A), 9 e [0, 27r), on the 

C-plane for the p listed in (|45p 


Proposition 14.41 extends to |a| = 1 as follows: 

Proposition 4.7. 


[exp(i(9Z^(2))l{.r<oo}] = exp(i(9z(2)) (^/3i(e*®)^ 


z{l)-z{2) 


and 


^i(e*®) = |E(i^o) [exp(f(9Z^(2))l{^<oo}]| ^ r. 
Proof. The dominated convergence theorem implies that 

(9 ^ E^ [exp(i(9Z^(2)l{^<oo}] 


(62) 

(63) 


(64) 
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is continuous in 9. Proposition 14.11 implies that for each fixed 0 the value of this function 
equals either 

exp(i 6 ' 2 :( 2 )) (65) 

or 

exp(i 6 'z( 2 )) (^/32 ^ ^ 

(l58l) and m imply that these two expressions are continuous as functions of 6 and they 
are never equal. Then (1641) must equal one or the other for all 9 and therefore if one can 
verify that (|62]) equals (l 6 ^ for a single 9 then the equality must hold for all 9; (|M|) asserts 
the desired equality at 0 = 1; (1521) follows. Set z = (1,0) in (|U2]) and take absolute values of 
both sides: 




lE(i,o) [exp(i6»Z^(2))l{^<oo}] 


I • I is convex; then Jensen’s inequality gives 

^ ®(1,0) [1{t<oo}] “ 

where the last equality is (IMI) with z{l) = 1 . □ 

For two tandem queues, the narrow shaded region of Figure 0] (let’s call it TZ) shows 
the set of parameter values that violate (|56|) (the shaded triangle containing 7Z shows the 
parameter values of stable tandem walks, i.e., the region X < fii, fi 2 where A = 1 — — ^ 2 )- 

The image of 0 —> A(e*^), 9 e [0, 27r), becomes a self intersecting curve on C when (1561) fails 



Figure 4: The set of parameter values violating (|56l) for the two dimensional tandem walk 
(the narrow shaded region on the left) 

and, e.g., the proof of Proposition 14.51 (or that of its adaptation to the parameter values in 
TZ) will require a study of )R(A) over TZ] we leave this analysis to future work. 
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4.4 log-linear harmonic fnnctions of Y 

Let us rewrite (1351) separately for the boundary 62 and the interior B° — d 2 - 

V{y) = J]V{y + v)p{v),yeB°-d 2 , (66) 

veV 

y{y) = V{y)lJ^ 2 + Yi V{y+ v)p{v),y e d 2 n B°. (67) 

vgV,v(2)j^-1 

Any g e Tz satisfies ([66]) (because (f66]l is the restriction of (|39|) to — ^ 2 ); (l66|) is linear 
and so any finite linear combination of members of Tz continues to satisfy (1661) . In the next 
two subsections we will show that appropriate linear combinations of members of Tz will 
also satisfy the boundary condition (1671) and define harmonic functions of Y. 

Parallel to the definitions in the previous section, we will define characteristic polynomials 
and surfaces for Y. The constrained process will have a pair of these, one set for the interior 
and one set for the boundary ^ 2 . The interior characteristic polynomial for Y is by definition 
that of Z, i.e., p and its interior characteristic surface is B. The characteristic polynomial 
of Y and its characteristic surface on the boundary 82 are defined below. 


4.4.1 A single term 

Remember that members of Tz are of the form [(/I, a), ■] : z —> 0 -^( 2 ) a^d (/3, a) e B] 

these define harmonic functions for Z and they therefore satisfy (I66p . The simplest approach 
of constructing a T-harmonic function is to look for [(/3, a), •] which satisfies (|35p . i.e., which 
satisfies (fM!) and (IS7D at the same time. Substituting [(,0,0!), •] in ([U7)) we see that it solves 
(1671) if and only if (/?, a) e B also satisfies 

P2(/0, o) - (da = 0 (68) 

where 

P2(/3, o) = /3a I ^ p{v)(3y^yy‘^'> + /X2 |. (69) 

y«GV,i;(2)7^-l / 

We will call ([68|) “the characteristic equation of Y on 02” and p 2 its characteristic poly¬ 
nomial on the same boundary. p2 can be expressed in terms of p as follows: 


P2 = p(/3,a) - /3a ^ - p2 

\veV,v{2) = -l 

= p(/?,a) -/5a (p{2,0)- +p{2,l)- - fi2^ . (70) 

\ a a / 

Define the boundary characteristic surface of Y for 82 as B 2 = {(/3,a) e : p2(/3,a) = 

0 }. 

Then, (/3, a) must lie on 77 n B 2 for [(/3, a), •] to be a harmonic function of Y. Suppose 
(/3, a) e B n B 2 , i.e., P2(/3, a) = p(/3, a) = /3a; and ,0, a ^ 0. Then, by (I7n|) 

0=p(2,l)-+p(2,0)^-/r2, 
a a 

/5 = ^^^(/^2a-p(2,1)). (71) 
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Substituting this back in (ilil) implies that a must solve 


P2(«) = 0> 


where 



■p(0,l) 


p(2,0) 

M2,0) 


/ _ , p(l,0M2,l) _ _ p(l, 2)p(2,1) 


P(2,0) 

P(1,0M2,1) 

M2,0) 


P(2,0) 


a 



(the superscript r stands for “reduced.”) Then [(/3,a), •] is a harmonic function of Y if and 
only if a is a root of P 2 and fi is defined by (1711) . The functions z —> 1 and z —>■ 0 are 
harmonic functions of Y of the form z ^ then two of the roots of P 2 are 0 

and 1 (that 0 is a root also directly follows from the form of P 2 ). It follows that the third 
root is 


p(0,l) 


n = 


P(2,l) 

P(2,0) 


1 + 


P(b0)p(2.1) 

P(2,0) 


1^2 


M2 / 

P(2,0) V 


p(i.o) 

^2 TO 


p(l,2 


This quantity is always less than 1 if the first queue is stable: 
Lemma 1. ui < fii if and only if ri < 1. 

Proof ri < 1 is equivalent to 


p(0,l) 


K2,l) 

p(2,0) 


(1 - ^ 2 ) < p(l,0) + 


2p(2,l)p(l,0) 

p(2,0) 


Ml, 2) 


M2,1)M1,2) 

p(2,0) 


substitute yii +p(0,1) +p(0, 2) for 1 — /U 2 , multiply both sides by p{ 2 , 0) and cancel out equal 
terms from both sides: 


P2P{0, 1) + P{2, l)p{0, 2) < mip(2, 0) + p{2, l)p(l, 0) (72) 

P2P{0, 1) + p(2, l)p(0, 2) < p{l, 2)p{2, 0) + p2P{l, 0) 

p(0,1) + 2) < Ml, 2)^^ + p(l, 0); 

P2 P2 

divide both sides by 1 — to get < pi. This establishes the “only if” part of the 

statement of the lemma. The last sequence of inequalities in reverse gives the “if” part. □ 

And thus we get our first nontrivial harmonic function for Y: 

Proposition 4.8. The function 

is a harmonic function of Y where 

/3(a) = Q) {P2a - M2,1)) 
is the right side of m- Furthermore 0 < /3(ri) < 1. 
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Proof. It remains only to prove the last part of the proposition’s statement. ri < 1 implies 
ri //2 < fJ -2 = p(2,0) + p{ 2 , 1) or, what is the same, 


this is /3(ri) < 1. The inequality f3{ri) > 0 turns out to be true for all p as long as p{2, 0) > 0 
and follows from a sequence of inequalities similar to ()72h . □ 

4.4.2 Two terms 

Define the boundary operator D2 acting on functions on and giving functions on 82'- 
D 2 V = g, V -.1? 

giy,0)=ip2+ 2 p(u)I7((y,0) + u) I - l/(y,0), y 6 Z; 

\ veV,v{ 2 )^-l ) 

(if V is defined on a subset of Z^ one may extend it trivially to all of Z^ to apply D 2 ). D 2 
is the difference between the left and the right sides of (1671) and gives how much V deviates 
from being y-harmonic along the boundary 82'- 

Lemma 2. D 2 V = 0 if and only ifV is a harmonic function ofY on 82 - 

The proof follows from the definitions involved. For (/3, a) e and /3,a 0 

[D2 ([(/3, a), •])] (y, 0) = (^^P2(/3, a) - 

where the left side denotes the value of the function D 2 ([(/3, a), •]) at (y, 0), y e Z. For 
(/3,a) e %, p(/3,a) = /3a; this, the last display and ([70]) imply 

[D 2 ([(/3, a), •])] (y, 0) = [p 2 - (pi2, 0)^ + p{2, 1)1) ) py (74) 

if (/3,a) e PL. One can write the function (y, 0) —> /3^ as [(/3, a),-Jlgj = [(/3,1),-Jlgj; in 
addition, define 

C(/3, a) = p2- (p(2, 0)1 + p(2, 1)1) . (75) 

\ a a/ 

With these, rewrite ([7i|) as 

D2([(/3,a),-]) =C(/3,a)[(/3,l),-]|a,. (76) 

The key observation here is this: D 2 ([(/3,a), •]) is a constant multiple of [(/3,1), ■]\d 2 - This 
and the linearity of D2 imply that for 

ai Y a 2 , (/3, ai), (/3, 02 ) e 77, (77) 

[(/3,ai),-] and [(/3,a2),-] can be linearly combined to cancel out each other’s value under 
D 2 . We will call (/3,ai) and (/3,a2) conjugate if they satisfy (i77|) . An example: the two end 
points of the dashed line in Figure [7] are conjugate to each other. 
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Because the characteristic equation (|44l) is quadratic in a, fixing /3 in (1441) and solving 
for a will give a conjugate pair (/3 ,q:i) and (/3 ,q: 2) satisfying 


/3 -p(l,0)/3^ -p(0,l) 
p(0,2) +/3p(l,2) 


(78) 


for most /3 e C; the next proposition uses these conjugate pairs and the above observation 
to define harmonic functions of Y: 


Proposition 4.9. Suppose that for ft e C, /3 ¥= 0, p{0, 2) + ftp{l, 2) ^ 0. Then 

hf} = C{P, a2)[(/3, oi), •] - C{P, ai)[(/3, 02 ), •] (79) 

is a harmonic function ofY. 

Proof By assumption (ft, oi), {ft, 0 : 2 ) are both on P and therefore [(/3, oi), •] and [(/3, 02 )) •] 
are harmonic functions of Z and in particular, they both satisfy (|66l) . Then their linear 
combination hjs also satisfies (f66]l . because (f66]l is linear in V. It remains to show that 
solves (1U7|) as well. /3 ^ 0 implies 01,02 ^ 0,1. Then (TTUl) implies 

D 2 {hp) = C{P, 02)T>2([(/3, ai), •]) “ C{(3, ai)D 2 {[P, 02 , •]) 

= C(/3,02)C(/3,oi)[(/3,l),-]|a, -C'(/3 ,oi)C'(/3,02)[(/3,1),-]|5. 

= 0 


and Lemma [2] implies that hjs satishes (|67p . □ 

Remark 4. If we set /3 = /3(ri) in the last proposition hjs reduces to a constant multiple of 

dZl. 

With Proposition 14.91 we define our basic class of harmonic functions of Y: 


Yy = {/i/3,/3 ^ 0,p(0, 2) + /3p(l,2) ^ 0}. 


(80) 


Members of Ty consist of linear combinations of log-linear functions; with a slight abuse 
of language, we will also refer to such functions as log-linear. 


Lemma 3. Suppose p(0,2) -I- /3p(I,2) ^ 0 so that (f78]l makes sense. Suppose further that 
(/3,ai) e P, (ft, 02 ) e P are conjugate with Q:i,a 2 ^ 0. Then (1781) is equivalent to 


p(2,0)/32+p(2,1)/3 
p(0,2)+/3p(I,2) 


(81) 


Proof. By ([78|) 


as-iiP - p(l, 0)/3^ - p(0,1)) _ 

a3-i(p(0,2) + /3p(l,2)) 

as-iP - p{l, 0)/3^ - p{0, 1) - g^O, 2) - /3a|_.p(l, 2) 
«3-i(p(0,2) -t Pp{1,2)) 


{a 3 -i,P) e P implies 

^ 1 p{2,0)P^+p{2,l)P 

as-i p{0,2) + Pp{1,2) 

□ 
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Define 


( 82 ) 


^ l p(2,0)/j2+p(2,l)/3 

a p(0,2) + /3K1,2) 


We can write m as 


Ui = a{a3-i). 

The map a is invertable (it is a multiple of a~^) and its inverse equals itself. Thus, conjugacy 
is symmetric: if (/?, 02 ) is conjugate to (/3,a2)j then (/3,ai) is conjugate to (/3,a2)- We will 
sometimes refer to cx as conjugator. 


4.5 Graph representation of log-linear harmonic functions of Y 

Figure [5] gives a graph representation of the harmonic functions developed in the last sub¬ 
section. Each node in this figure represents a member of J-z- The edges represent the 


2 



Figure 5: The harmonic functions of Y 


boundary conditions; in this case there is only one, (I67|) of ^ 2 , and the edge label “2” refers 
to ^ 2 . A self connected vertex represents a member of J-z that also satisfies the 82 boundary 
condition (f671) . i.e., z —> of Proposition ITHl the graph on the left represents 

exactly this function. The “2” labeled edge on the right represents the conjugacy relation 
()81h between ai and 02 , which allows these functions to be linearly combined to satisfy the 
harmonicity condition of Y on 82 ■ 

4.6 5i?-determined harmonic functions of Y 

Our task now is to distinguish a collection of 5i?-determined members of J-y- This collection 
will form a basis of harmonic functions with which we will approximate/ represent the rest 
of the 5i?-determined functions of Y. 

Proposition 4.10. Let ai, 02 0^4 /3 be as in Proposition |/.g[ If 

|/3| < 1, |q;i|, |q;2| ^ 1 (83) 

then of (17^ is dB-determined. 

Proof. By Proposition 14.91 Hr is a T-harmonic function; (I83p and its dehnition (|79l) imply 
that hg is also bounded on B°. Then = h^(TT-AC„AA:) is a bounded martingale. This, 
Proposition 12.21 and the optional sampling theorem imply 

hgiu) = lEy [/i^(l(-)l{r<Cn}] ”1" \_^dO^Cn)^ {CnsSr}] ) y ^ BL (84) 

y/^(l) = n for T > Cn- This and (f83]l imply 

Jin^E, ^ ^Im /?" = 0. 
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This, lim„ = oo and letting n ^ oo in (|84p give 


hpiy) = Ej; [hp{Yr)l {r<oo}] ; 


i.e, hfi is 5i?-determined. □ 

Remark [H Proposition 14.81 and the last proposition imply 

Proposition 4.11. The harmonic function (|73p is dB-determined. 

What Proposition 14.101 does is it gives us a collection of basis functions for which the 
Balayage operator T(^b°Y is extremely simple to compute; these functions play the same role 
for the current problem as the one which exponential functions do in the solution of linear 
ordinary differential equations or the trigonometric functions in the solution of the heat and 
the Laplace equations. Let us rewrite Proposition I4.1UI more explicitly. Suppose ai, 02 and 
(3 are as in Proposition 14. lOl Define 

f{y) = a 2 )af^ - Cif3, ai)af^ ,yeZ^ 

Then, Proposition 14.101 savs 

E, [/(W)l{r<oo}] = (c{/3,a2)af^ - C{f3,a,)af^) 

= /3''(i)-^(2)/(y), yeB. 

Proposition 14.101 rests on the condition (I83h : Proposition 14.131 below identifies a set of con¬ 
jugate pairs (/3,ai) and (/3,a2) on TL satisfying ([831) . 


4.7 A modified Fourier basis for T{b°y 

Let’s go back for a moment to the problem of evaluating the Balayage operator of the 
unconstrained process Z for the set {B^Y. For a bounded function / on dBz, this is the 
operator mapping / to the harmonic function 

[fiZr)l{r< 03 }] ,ze Bz] 

in the present subsection we will write T for T(b'^Y Fourier basis functions 

fa ; dBz C, /«(y, y) = y G Z, |a| = 1, (85) 


we already know how to compute T{fa) (given by ([62])) and T is linear. One can use these 
to evaluate T more generally in three related ways. First, Fourier series theory tells us that 
if / is I 2 , i.e., if Xiysz IfiYiY)]'^ < 00 , it can be written in terms of the Fourier basis functions 
thus: 


/(y,y) 


1 



00 

mYy^deJ{0)= 2 /(y,y)e-*^' 

y=-co 


( 86 ) 


Fubini’s theorem now implies T{f) = ^ /(^)7”(/e io^dO^ and one can construct approxi¬ 

mating sequences by truncating the sum in (I86p . 

Second, when interpreted as a function of a, the formula (|62p gives the characteristic 
function of the distribution T. Its inversion would give the distribution T itself. 
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Third, we can first replace / with its periodic approximation defined as follows 


F{y,y) 


/(y,y), if |y| 

/(y mod N,y mod N), if |y| ^ N. 


As N increases, |T(/^) — T(/)| will converge to 0. Because it is periodic, has a unique 

, 2A^ I 1 i( ^ 

Fourier representation of the form fP = Xlfc=o where Ok = e \2 n+i)_ Then T{f^) = 


How should one proceed to build a parallel theory for the constrained process Y1 The first 
obstacle to the above development in the case of Y is that the Balayage of the Fourier basis 
functions is not simple to compute, i.e., we don’t know a simple way to compute T(^B°)‘={fa)- 
But Proposition 14.101 savs that if 


a = Q:(/3i(a), a)| < 1, C{f3i{a),a) ^ 0 


(87) 


then the Balayage of the perturbed Fourier basis function 

/ \ y C^ (Yj , /.y 

(y,y) ^ - {a y,ye Z+, 

o (pr (ck), cy ) 


is simple to compute and is given by 

y 


ft(„),(.)-,(2) (,,( 2 ) _ ,yeB. 

V Y{pi[a),a ) 


( 88 ) 


(89) 


One can interpret (l 88 l) as a perturbation of the restriction of (l 8 ^ to dB, because \a'\ < 1 
implies that these two functions equal each other for y large. Below in Proposition 14.131 we 
identify a set of a’s for which (I87p holds and, therefore, for which the image of () 88 h under the 
Balayage operator is given by ()89p . This will require further assumptions on p; in particular, 
we will assume (IHTl) for a = 1 : 

a{r, 1) < 1, (90) 

C'(r,Q(r, 1)) ^ 0, (91) 

where, as before, r is the input/output ratio (fT71) . 

Comments on these assumptions: 

1. If p(2,1) is large enough, a{r, 1) can exceed 1 even when the stability assumption (jl4p 
holds. For such networks, we cannot check whether hr is dB-determined by an appli¬ 
cation of Proposition [4.101 Furthermore, even if hr is dB-determined, [(r, Q:(r, 1)), •] 
will dominate hr for large y{2) and one can no longer think of hr/C{r,cx{r,l)) as a 
perturbation of the function 1 on dB. 


2. The condition (|9ip implies that the log-linear function [(r, cx{r, 1)), •] is not T-harmonic. 
For two dimensional tandem walk, reduces to pi ¥= ^ 2 - One can treat the case 
hi = h 2 by writing it as the limit hi h 2 - This limiting procedure introduces 
harmonic functions with polynomial terms, see subsection 19.41 

With pnp and (f^ we are able to take a = 1 in (f88l) . Subsection 14.4.11 implies that for 
a e and a ^ 1, C(/3i(a), a) ^ 0. Thus, for such a only the first condition in ([87P is 
nontrivial. In the next proposition we will show that (19011 implies |q;(/3i(q;), a)| < 1 for all 
|a| = 1. To simplify its proof, we will further assume p(0, 2) = 0. Treating p(0, 2) 7 ^ 0 seems 
to require a more refined analysis of /3i as a function of p and a, a task we defer to future 
work. 
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Proposition 4.12. Suppose p{0,2) = 0, (IMt) and (l52T) hold; for p{l, 2) = 0 the system 
becomes trivial, so we will also assume p{l, 2) ^ 0. Then 


\a{Pi{a),a)\ < 1 (92) 

for all |a| = 1. 

Proof. The definition of ct and p(0,2) = 0 imply cx{P,a) = c{a)/a where 

c(a) = (p(2,0)/3i(a)+p(2, l))/p(l,2). 


Then 

|Q:(/3i(a), a)| = |c(a)| (93) 

because |a| = 1. As a varies on c(a) defines a closed curve in C that is symmetric around 
the real axis. ()63p implies that this curves is contained in a circle centered at the point 
p(2, l)/p(l, 2) with radius rp(2,0)/p(l, 2), which in turn is contained in the circle centered 
at the origin and with radius (p(2,1) + rp(2, 0))/p(l, 2) = a.{r, 1) which by assumption (IMl) 
is less than 1; then |c(q;)| < 1. This and P5]) imply (1^21) . □ 

We have assumed p(0, 2) = 0 only to simplify the above proof; Figure [6] shows an example 
with p(2,0) ^ 0 where again |q:| < 1. 



Figure 6: The graph of cx{fd, a), a = exp(z0) on the C-plane as 6 varies in [0, 27r) for the p 
of (li5|) 

Remark 5. If p(2, 0)/p(l, 0) = 1, a = —1 violates (15211 and the last proposition is not 
applicable at a = —1, because /3i is not well defined there. But in that case the single 
root of the affine (14611 will take the place of fi\ above. With this modification, the above 
argument works verbatim for p(2,0)/p(l, 0) = 1 as when p(2,0)/p(l, 0) 1 and from here 

on we assume that a similar modification is made when a violation of (I52p occurs. 

Proposition 4.13. Assume (fTT|l . ([90]), (IMIl and (f56]l . Then there exists 0 < R < 1 such 
that for all a e C with R < |a| ^ 1, is a dB-determined harmonic function ofY. 
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Proof. Proposition 14.101 and Proposition 14.12l imply that is 5i?-determined harmonic 

function of Y for a e S^. a and /3i are continuous functions. This, compactness of and 
([63]) imply 

|/3i(a)|, \oL{l3i{a),a)\ < 1 

for |a| e [R, 1] where R is sufficiently close to 1. This and Proposition 14.101 now imply the 
statement of the proposition. □ 

5 Analysis of y, d> 2 

Now we would like to extend some of the ideas of the previous section to d dimensions. Our 
path will be this: each of the graphs shown in Figure 0 and its corresponding equations 
define a harmonic function of Y. We will develop the graph representation in d dimensions 
and show that any solution of the equations represented by a certain class of graphs defines 
a harmonic function of Y. 

For a e we will index the components of the vector a with the set M = A/q — {0,1}, 
i.e., a = (Q;(2),a(3), ...,a{d)). The members of M are exactly the constrained coordinates of 
Y. For a c AA the constraining map vr sets the following set of increments of T to 0 on da- 

Va = {v. V eV, v[j) Y -1, for j e a}. 

Rewrite (I35|) as 

y{y) = + Yi ^ danO,a^J\f. (94) 

jea veVa 

Set 

[(/3, a),z\= / 3 Ai)-ye^Ai) (95) 

jej/ 

[(/3, a),z] is log-linear in z, i.e., log([(/3, a), z]) is linear in our goal is to construct Y- 
harmonic functions out of linear combinations of these functions. 

Define the characteristic polynomial 


Pa(^,a; 

1 = ( 2 

+ 2/^i I 

(96) 


\veVa 

jea } 


the characteristic equation 





Pa(^,a) = 

1, 

(97) 

and the characteristic surface 




Ua 

= {(/3,a)eC'':p 

aid,a) = 1 } 

(98) 


of the boundary da, a cz M. We will write p and PL instead of p 0 and 'H 0 . 
Generalize (HQI) to the current setup as 


Pa(^,a) = p(/3,a) - I 2 p(u)[(/3,a),u] - . (99) 

\veV^ jea J 

Pa is not a polynomial but a rational function; to make it a polynomial one must multiply 
it by dWjeN what we did when we gave the two dimensional versions of these 
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definitions in ([jH]) and (l6^ . For d > 2, the /SfljeA^aO) multiplier complicates notation; 
for this reason we omit it but continue to refer to the rational (j96l) as the “characteristic 
polynomial.” 

The argument in subsection 14.4.11 continues to work verbatim for general d and gives 

Lemma 4. Suppose {/3,a) e T-L n 'Hi. Then [(/3,a),-] is Z-harmonic (i.e., Y-harmonic on 
fly ) and Y-harmonic on di. 

Our next step is to extend the content of subsection 14.4.21 to the current setup. Begin 
with the operator Da acting on functions on and giving functions on da- 

DaV = g, 

aiz) = ( 2 ^“+ S pivWi^ + v)] -v{z). 

\jea veVa } 

Lemma 5. DaV = 0 if and only ifV is Y-harmonic on da- 

The proof follows from the definitions. Next generalize C of (17511 to 

C{j,/3,a) = Hj - ^ p{v)[{l3,a),v]. (100) 

VGV,v{j)= — l 

For a e and a ^ M define a{a} e as follows: a{a}|ac = a\a<^ and a{a}\a = 1; 
If a = {z}, we will write a{i}, instead of a{{z}}. For example, for J\f = {2,3,4}, a = {4} 
and a = (0.2, 0.3,0.4), a{a} = (0.2,0.3,1). We will use this notation in the next paragraph, 
where we define the conjugacy of points on 7^ in d dimensions and in the next section where 
we apply the results of the present section to d-tandem queues. 

For i e M, fix and multiply both sides of the characteristic equation (ITZll by 

a;(z); this gives a second order polynomial equation in a{i). If /3 and are such that 

the discriminant of this polynomial is nonzero, we get two distinct points (/3, ai) and {l3, 02 ) 
on Ti which satisfy 


“iIn-w 

Q;i(i) + a2{i) 


I]^(j)=i[(/3,Q;i{z}),z;] ■ 


( 101 ) 

( 102 ) 


The sum in the numerator on the right side of (11021) is over v such that v{i) = 0; this and 
(jlOip imply that (110211 remains the same if we replace [(/ 3 ,q;i),u] in the numerator with 
[{P,a 2 ),v] or [{fd,a 2 {i}),v] = [{(3,ai{i}),v]. Rewrite (I102p as 


a2{i) 


1 - 2 ^ 0=0 

i;„(i)=i[(/3,ai{i}),u] ai(z) i;„(i)=i[(/3,aijf}),u] 


(103) 
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Now keep 02 ( 1 ) on the left and repeat the same computation to get 


1 Ilv{i)=-i[{l3,a2{i}),v] 

«2(i) 'Lv(i)=i[il3,a2{i}),v] ■ 


( 104 ) 


We will call (/3,ai) ¥= (13,02) £ 33 i-conjugate if they satisfy (|101l) and any of the equivalent 
()102p . (llOOh and (1104^ . Based on these, generalize the conjugator a as cx(i,(l3,ai)) = 02 
where 02 is defined by (llOip and (jl03p . 

Our next proposition generalizes Proposition 14.91 to the current setup. In its proof the, 
following decomposition, which (I99p implies, will be useful: 


Da([(l3, a, ODCz) = I XI “ X 'P(v)[(l3, a), u] 1 [((3, a), z] (105) 

\jGa veVS j 

= X(^1“ X P{v)[{l3^a),v]\[{l3,a),z] 

jea y veV,v(j)=—l j 

.X;0)([(/3.a.-)])W (106) 

jea 

for z e da and (/3, a) e T-L. 

Proposition 5.1. Suppose that (f3,ai) and (( 3 , 02 ) are i-conjugate and C(i, (3,aj), j = 1,2 
are well defined. Then 

hp = C(i, 13, a2)[il3, oi), •] - C{i, fi, ai)[(/3, 02 ), •] 

is Y-harmonic on di. 

Proof. The definition piOOp of C, pi05p and linearity of Di imply 

A(V) = C{i,fi,a2)C{i,fi,ai)[{fi,ai),-] - C{i, fi,a2)C{i, fi,ai)[{fi,a2), ■] 


(llOip implies \{l3,ai),z) = \{f 3 ,a 2 ),z) for z e di and therefore the last line reduces to 

= 0 . 

Lemma [5] now implies that hp is T-harmonic on di. □ 

To generalize the graph representation of the previous section we will need graphs with 
labeled edges; let us denote any graph by its adjacency matrix G. Let Vc, a hnite set, denote 
the set of vertices of G . Each edge of G will have a label taking values in a finite set L. For 
two vertices i Y j, G(i,j) = 0 if they are disconnected, and G(i,j) = Hf an edge with label 
I e L connects them; such an edge will be called an Ledge. As usual, an edge from a vertex 
to itself is called a loop. For a vertex j e Vg, G(j,j) is the set of the labels of the loops on 
j. Thus G(j,j) c L is set valued. 

Definition 5.1. Let G and L he as above. If each vertex j e Vg has a unique l-edge (perhaps 
an l-loop) for all I e L we will call G edge-complete with respect to L. 
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We say G is edge-complete with respect to Y if it is so with respect to J\f = Mq — {0,1} 
(remember that J\f is the set of constrained coordinates of Y). If we jnst say “edge-complete” 
we mean “edge-complete with respect to y.” 


Definition 5.2. A Y-harmonic system consists of an edge-complete graph G with respect to 
N, the variables {/3,aj) e Cj eC, j e Vg, and these equations/constraints: 

1 . (/3, Oj) G Cj e C — {0}, j e Vg, 

2 . Ui Y aj, ifi Y j,ij e Vq, 

3. ai,aj are G{i,j)-conjugate ifG{i,j) YO, i Y j,i,j e Vq, 


4- 

5. 


^ G{G[i,j),l3,aj) 

(/3, a/ e Ui for all I e G{j,j),j e Vq. 


ifG{i,j) Y 0, 


(107) 


Proposition 5.2. Suppose that a Y-harmonic system for and edge-complete G has a solu¬ 
tion; then 


hG= Yi 


(108) 


jeVa 


is a harmonic function ofY. 


Proof. All summands of ho are Z-harmonic and therefore y-harmonic on fly because (/3, a/, 
j e Vg, are all on the characteristic surface PL. It remains to show that Pg is y-harmonic 
on all da n Hy, a M and a V 0. We will do this by induction on |a|. Let us start with 
|a| = 1, i.e., a = {1} for some I e M. Take any vertex i e Vg; if I e G{i,i) then {(3, a/ e PLi 
and by Lemma H] [(/3,aj),-] is y-harmonic on di. Otherwise, the definition of a harmonic 
system implies that there exists a unique vertex j of G such that G{i,j) = 1. This implies, 
by definition, that (/3, Oi) and {(3, a/ are /-conjugate and by Proposition 15.11 and (I107p 


Cj[(/?,«*),•] -t Cj[(/3,aj), •] 

is y-harmonic on d;. Thus, all summands of Pg are either T-harmonic on di or form pairs 
which are so; this implies that the sum Pg is y-harmonic on di. 

Now assume Pg is y-harmonic for all a' with \a'\ = k — 1; fixai^A/” such that \a\ = k 
and i e a; by (11061) 

Da{hG) = Da-{i}{hG) + D/hG). 

The induction assumption and Lemma [5] imply that the first term on the right is zero; 
the same lemma and the previous paragraph imply the same for the second term. Then 
Dai^G) = 0; this and Lemma 0 finish the proof of the induction step. □ 

Proposition 5.3. Let (/3,aj), Cj, j e Vg, he the solutions of a Y-harmonic system for an 
edge-complete G and let hp he defined as in (llU8p . If 

|/3| < 1, |aj(/)| ^ 1, ieM, 


then Pg is dB-determined. 

The proof is identical to that of Proposition 14.101 
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5.1 Simple extensions 

One can build, from the solution of a given harmonic system for Y, solutions to related 
harmonic systems for higher dimensional walks which are, in some natural sense, extensions 
of Y. The construction will depend on what we mean by an “extension.” One possibility is 
that of a simple extension whose definition follows. 

So far, we have taken p to be a (d + 1) x (d + 1) real matrix, i.e., p e where 

d + 1 = |A/o|. To dehne “simple extensions” it is more convenient to take the set of nodes A/q 
to be an arbitrary set with d + 1 elements containing 0 and p e A/+ is, as before, 

1 A/*^xA/*^ 

A/q ~ {0}- Suppose Aq A/q and pi e M_,_° ° is another matrix of jump probabilities. 

Define p' e as follows 

P'ihj) = PiiiJ) (109) 


if i e Afo, j e Af+, 

p{i,0) = pi{i,0) + Yi Pi(bi)- 
jeA/g-A/b 

Definition 5.3. We say that pi is a simple extension of p if 

P' = ( 2 ) P^P' ^ 

\i,jeAfo / 

Pi{hj) =0 i/i e A/'o - A/j) and j e M+. 
An example: take Aq = {0,1,2}, A/g = {0,1, 2, 3,4}, and 
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Figure [7] shows the topologies of the networks corresponding to p and pi. 


( 110 ) 


( 111 ) 

( 112 ) 



Figure 7: Two networks, second is a simple extension of the first 
Next define the “edge-complete extension” of a given edge-complete graph: 
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Definition 5.4. Let G be an edge-complete graph with respect to a finite label set L . Its 
edge-complete extension Gi with respect to another set of nodes Li ^ L is defined as follows: 
Vgi = Vg and 

Gi{i,j) = G{i,j),i j,i,j e Vg 
G i(j,j) = G(j,j) u (Li - L),j e Vg- 

To get Gi from G one adds to each vertex of G an /-loop for each I e Li — L. Then if 
G is edge-complete with respect to L, so must be Gi with respect to Li. Figure [8] gives an 
example. 




Figure 8: An edge-complete graph with respect to {2} and its edge-complete extension to 
{2, 3,4, 5} 


Suppose is a simple extension of Y ; the next proposition explains how one can con¬ 
struct harmonic systems (and their solutions) for Y^ from those of Y. 

Proposition 5.4. Let Y^ be another constrained process defined using the construction (I22p 
(in particular i = 1 and only the first coordinate ofY^ is unconstrained and its remaining 
coordinates are constrained) and such that its matrix of jump probabilities pi is a simple 

extension of the p matrix ofY. Let Gi and Gq be edge-complete graphs for Y^ and Y such 

that Gi is an edge-complete extension of Gq. Suppose (/3, ak), c^, k e Vgq, solve the harmonic 
system associated with Gq. For k e Vg^ = Vgq define a\ as follows 

alW = ak, (113) 

ailAfi-AT = fi, (114) 

where ^ N is the set of constrained coordinates of the process Y^. Then {(5,a\),Ck,k e 
Vgi ; solve the harmonic system defined by Gi. 

The definition (111411 assigns the value j3 to the new components of a\ coming from the 
new dimensions of the simple extension; this corresponds to ignoring the new dimensions 
when we compute the log-linear function \ffi,a\), •] at the increments of Y^ (see (111611 and 
()117p below). The following proof lays down the details of this observation. 

Proof Set A/)| = u {0,1} and A/":} = JY^ u {!}. By assumption, {fi,ak), Ck, k e Vgq, 
satisfy the five conditions listed under Definition 15.21 for G = Gq. We want to show that this 
implies that the same holds for (/3,a^), Ck, k e for G = Gi. Let Vq denote the set of 
increments of Y^, ej, j e JY)) the unit functions on A/":} and let ej be the 0 function on the 
same set. (jll2p implies that we can partition Vq as follows: 

VqI = V^ u Vs^ u Vl, (115) 

4^1 = {e\ + e],-e} -e},-e- -he],zGA/'u {0},j eM}, 

V 2 = {e\ + + e],ie JY,j e [jYi - JY) u {0}}, 

Yq = {-e{ + e],i,j e {JYi - JY) u {0}}. 
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Parallel to this is the following partition of Vq: 


Vo = Vi u V2, 

Vi = {ei + ej, -Ci - ei, -e* + ej,iej\f u {0}, j e Af}, 

V 2 = {ei, -ei,ie Af}. 

Fix any k e (I113p and (11141) imply 

[{I3,al),v^] = [{l3,ak),v\ (116) 

for all e u V 2 and where v = (11141) implies 

= 1 (117) 

for e V 3 . Let denote the characteristic polynomial of and let 71^ denote its 
characteristic surface; we would like to show {j3,a\) e 77}, i.e., p^(/3,a^) = 1. By (11121) and 

(HEI) 

pH/ 3,afc)= Ya Pi(^H[(/3,afc),'f^H + Y Pi(^H[(/3,(US) 

v^eVl v^eV2 

+ Y Pi(^H[(/3,afc),^^H- 

IjlsVg 

(I109P and (11161) imply 


= 2 H('^)[(/3,afe),^^] + Y Pi(^H[(/3,afc)>^^H + Y (119) 

^gVi ,;1sV| YeVl 

For any 6 V 2 , (I114p implies 

[(/3, Oil),v^] = [(/?, ak),v\ (120) 

where v = this implies that the second sum on the right side of (|119p equals 

Y P (^)[(/3,«fc),^]- 

veV 2 

Substitute this back in (|119l) to get 

pH/3,al) = 2H(^)[(/3,afc),^^] + Y 

veV uisVj 

which, by (lllip . equals 

= (2p'(^)) + Y 

\veV ) veV 


(/3, ak) G 77, piOOp and pilOp now give 

= Y + Y = 1, 

v^eVlvjVl v^eVl 
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i.e., indeed, e T-L^. 

Let us now show that the third part of the same definition is also satisfied. Fix any i ¥= j 
with = I e M. We want to show that {/3,aj) and (/3,aj) are /-conjugate, i.e., that 

they satisfy (jlUip and (|lU3p : 

( 121 ) 


= 7J7 


1 


( 122 ) 


l>vm)=iPi{v)[{P,a}{l}),v^] ■ 

By dehnition = I when G{i,j) = /; and this implies, again by definition, that Oj 

and aj satisfy (jlDip and (110311 . (112111 follows from (110111 for Oi and Uj, (I113p and (111411 . For 
I e Af aj{l) = aj{l) and al(/) = «*(/). Then to prove (|122p it suffices to prove 


Pi{v)[{P, a,^/}), p{v)[{/3, ai{l}),v] 


(i)- 


Pi{v)[{P, aj{l}),v^] Y,vHi)=i Piiv)[{P, ai{^}),v] 


(123) 


This follows from a decomposition parallel to the one given for the previous part: let us first 
apply it to the numerator. 


2 Piiv)[iP,al{l}),v^] 

Yi Pi{v)[{P,al{l}),v^] 


Y Pi{v)[{P,a\{l}),v^] 

'ui(Z)=-l,'uisV] 


(IllOp and (I120p imply 


Y P\'^)[i/3,o:i{l}),v] + Y Piiv)[{P,ai{l}),v] 

v(l)=—l,veVi v(l) = — l,veV 2 

V^’gV / v{l)=—l 

A parallel argument for the denominator gives 


(124) 


•ui(0 = l / v{l) = l 


Dividing (|124ll by the last equality gives (I123p . 

The proof that the remaining parts of the dehnition holds for G = Gi is parallel to the 
arguments just given and is omitted. □ 


6 Harmonic Systems for Tandem Queues 

Throughout this section we will denote the dimension of the system with d; the arguments 
below for d dimensions require the consideration of all walks with dimension d ^ d. 

We will now dehne a specihc sequence of edge-complete graphs for tandem walks and 
construct a particular solution to the harmonic system dehned by these graphs. These 
particular solutions will give us an exact formula for Py{T < oo) in terms of the superposition 
of a hnite number of log-linear T-harmonic functions. 
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We will assume 


IJ-i ¥= fJ-jii ¥= j- (125) 

This is the analog of (|9ip for the d dimensional tandem walk. One can treat parameter 
values which violate (11251) by taking limits of the results of the present section, we give 
several examples in snbsection 19.41 

The p matrix of the tandem walk is as given in (jl5p . Then its characteristic polynomials 
will be of the form 


p{f3,a) = + l^ia(2) + ^ 

Pi{f3,a) = X-+ fiia{2) + fii+ \ pj 

p . aij) 


(126) 


where by convention a{d + 1 ) = /3 (this convention will be used throughout this section, and 
in particular, in Lemma [ 6 l (11271) and (11281) 1. The formula ()126l) for p implies 

Lemma 6. (/3 ,q:) eT-Ln^T-Lj {/3,o:) e V., pj = Pj e %, a{j + 1) = 

«(j); j 


The conjugators for the d-tandem walk are: 


a(^,(/3,«))(/) 
ol{1, (/3,q;))|^_p} 


1 a{3)iJ.2 

a( 2 ) 111 ’ 

1 a{l-l)a{l+l)iii 

. «(b W-l 




1 = 2 , 

2 < I ^ d, 


(127) 


For tandem walks, the functions C[j,j3,a) of (llOOl) reduce to 

C{j,P,a)=p,-p,^^^-^. (128) 

«(j) 

We define the edge-complete graphs d e {1, 2, 3,..., dj: 

Vg, = {au {d}, a cz {1,2, 3,..., d - 1}}; (129) 

for j e {avj {d}) n M define by 

Gd{a u {d}, a u {d} u {j - 1}) = j ii j - I i a (130) 

and 

Gd{a u {d}, a u {d}) = Af — a u {d}; (131) 

these and its symmetry determine Gd completely. Figure [9] shows the graph G 4 for d = 4. 
The next proposition follows directly from the above definition: 


Proposition 6.1. One can represent Gd+i as a disjoint union of the graphs Gk, k = 1,2 ,.., d, 
and the vertex {d -I- 1} as follows: for a cz {1, 2,3,..., A: — 1} map the vertex a u {k} of Gk 
to vertex a u {k,d + 1} of Gd+i- This maps Gk to the subgraph of Gd+i consisting of 
the vertices {a,k,d -I- 1}, a c {1,2, 3,..., A; — 1}. The same map preserves the edge structure 
of Gk as well except for the d + 1-loops. These loops on Gk are broken and are mapped to 
d -I- 1 -edges between and G'^_^_^. 
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Figure 9: G 4 for d = 4 


Figure [9] shows an example of the decomposition described in Proposition 16.II 
Define 


|a|-l a(j + l) 


< - n n 

<{i) = ^ 


/.i/ - A 


j = l l=a{j) + l ^“(1) 

1 if Z ^ a(l) 

Pa{j), if a(j) < ^ ^ a(j + 1 ), 
Pa{\a\) if I > a(|a|). 

f^a ~ Pa{\a \). 


(132) 


(133) 


I e J\f (remember that we assume that the elements of sets are written in increasing order; 
a(|a|) then denotes the largest element in the set). Several examples with d = 8: 


i5} 


= 1 . = ( 1 . 1 . 1 . l.PS./^S.Ps), 


-{ 3 , 6 } 


-{ 3 , 5 , 7 } 


= (-l) 


P5} 

/i4 — A /xs — A /rg — A , , 

- 5 “{3 6 } = {^A,P3,P3,P3,Pe,P6), 

PA - P3P5- P3 Pa - P3 ^ ^ 

2 Pi-\ M5 - A p&- \ P7 - \ 


PA - P3P5- P3 Pa - Pb Pi- P5 
“{3,5,7} = (1) 1.P3,P3,P5,/55,P7), 

6^{3| (f.l? P3 , P3^ P3^ P3 , P 3 ) .^^{s} (f.l.f.f.f.l.f)? 


(134) 


remember that we index the components of a* with A/"; therefore, e.g., the first 1 on the 
right side of the last line is Q!*gj( 2 ). 
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It follows from (jl32p and (11331) that 


c 


* 

a'u{di ,d2} 




n 


l=di-\-l 


idi- \ 
1^1 


a 


* 

au{di} 


= a 


* 

au{di ,d} 


for any a(|a|) < di < d2 ^ M and a cz J\f. These and Proposition [ 6 T] imply 
Proposition 6.2. For d < d and y e dB 

-( n E E <Wl,0,V^ (135) 

\l=d+l aeVa^ 

^ d 

Proposition 6.3. For d ^ d, let Gd he as in (I129P and (11301) . 
a CZ {1,2,3, ...,d — 1}, defined in (|132l) . solve the harmonic system defined by Gd- 


Proof. The first d components of a tandem walk is a simple extension of the tandem walk 
consisting of its first d — 1 components. This and Proposition 15.41 imply that it suffices to 
prove the current proposition only for d = d. 

Let us begin by showing ^pd, ® ■^+ characteristic surface Td 

of the tandem walk. We will write a* instead of a* r , the set a will be clear from the 
context. 

Let us first consider the case when a(l) > 1, i.e., when 1 a; the opposite case is 
treated similarly and is left to the reader. Then a*{l) = 1 for 2^1^ o(l)- By definition 
a*{i) = a*{i + 1 ) if a(j) < i < a{j + 1 ); these and = Pd give 


v{Pd,a*) 


Pd+ 2 Pj + Pa{l)Pa{l) + 2 Pj 

i=l jg(aii-{l...a(l)-l}) 


+ E 

jE(a-{a(l)}) 


a*{j + l) 
a*{j) 


+ Pd 


Pd 

a*{d) 


(where a^ = (A/+ —{d})—a) and in the last expression we have used the convention a*{d+l) = 
/3*; by definition (11331) a*(a(j + 1)) = Pa{j)^ 0 (*(o(j)) = Pa{j-i) therefore 


a(l)-l |a| 

= Pd + ^ Pj + X + E PJ + 'Ei P<i )~— - -^ ^“(l“l) 

j = l js(aii-{l-a(l)-l}) i=2 P«0-1) 

Pa(j)Pa(j)/Pa(j—1) Pa(j — 1) implies 

a(l)-l |a| 

= hd + 2 + A + ^ Pj + '^ Pa{j-l) + Pa{\a\) 

j = l js(a'i-{l-"a(l)-l}) j=2 

= Pd + E Pj ^ 2 Pj + Ep^ ^ 

j = l js(a=-{l...a(l)-l}) jea 
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i.e., {pd,a*) e U. 

If ai ¥= 02 take any z e oi — 02 (relabel the sets if necessary so that oi — 02 0). Let 

j be the index of i in oi, i.e., ai(j) = i. Then by definition, + 1) = Pi\ but i ^ 02 

and ()125p imply that no component of equals pi, and therefore 

This shows that a N satisfy the second part of Definition 15.21 

Fix a vertex a u {d} of Gd- By definition, for each of its elements I, this vertex is 
connected to o u {d} u {Z — 1} if / — 1 ^ o or or to o u {d} — {/ — 1} if / — 1 e o. Then 

to show that the a M satisfies the third part of Definition 15.21 it suffices to prove 

that for each a M and each I e avj [d] such that / — 1 ^ a u {d}, and 

are /-conjugate (remember that the graphs of harmonic systems are symmetric). For ease of 
notation let us denote a u {/ — 1 } by oi, by a*, by a\ and /3*^ = /3* by /?* 

(becanse we have assumed d = d, (3* is in fact equal to pd)- We want to show that {I3*,a*) 

and (/3*,a*) are /-conjugate. Let us assume 2 < I < d, the cases I = 2, d are treated almost 

the same way and are left to the reader. By assumption / e a* but I — 1 ^ a*. If / is the 

element of a, i.e., / = a(j); then a{k) = ai{k) for k < j, oi(j) = 1 — 1, a{k — 1 ) = ai{k) for 

k > j. This and (11331) imply 

(136) 

i.e., a* and a\ satisfy (llOip (for example, for d = 8 , 0*3 gj is given in (11341) : on the other hand 

®{3 5 g} (f 5 f) P3) P3) P5) P 6 ) P 6 ) Pg) tTrid indeed gjlAt—{ 6 } ®{ 3 , 5 , 6 }l.^~{ 6 } Definition 

(I133p also implies 

(/) Pa\{j) Pl—li (^ T 1) Paj(j-)-i) pi- (137) 

On the other hand, again by (jl33p . and by / — 1 ^ a, we have 


a*(/) = a*(/ — 1 ) = Pa[j-i) and a*{I - 1 - 1 ) = pi- 


Then 

1 a*{l-l)a*{l + l)pi 

*//\ Pl—\ 

a*{l) pi^i 

and, by (I137p this equals a*(/); thus we have seen a{l,{/3* ,a*)){l) = a*{l), where the 
conjugator cx is defined as in (I127p . This and (11361) mean that a\ = a.{l, (13*, a*)), i.e., 
{(3*, a*) and {(3*, a) are /-conjugate. 

Now we will prove that the a cz M, defined in (jl32p satisfy the fourth part of 

Definition 15.21 The structure of Gd implies that it suffices to check that 


_ C{l,pd,a*aJ 
<1 C{l,pd,cxt) 


(138) 


holds for any / e a such that / — 1 ^ o and ai = au{/ — 1}. There are three cases to consider: 
I = 2, I = d and 2 < / < d; we will only treat the last, the other cases can be treated 
similarly and are left to the reader. For 2 < / < d one needs to further consider the cases 
a(l) = / and a(l) < /. For b cz J\f, of (11321) is the product of a parity term and a 

running product of d — 6 ( 1 ) ratios of the form {pi — X)/{pi — Pa{j))- The ratio of the parity 
terms of o and ai is —1 because ai has one additional term. If a(l) = / then ai(l) = 1 — 1 
the only difference between the running products in the definitions of c* and c* is that the 
latter has an additional initial term {pi — X)/{pi — Pi~i) and therefore 


c; ^ Pi - Pl-l 
Pi-^ 
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Because I > 2 and 1 — 1^2, (ll.S.Sp implies a*{l) = 1, a*{I + 1) = pi, al{l) = pi-i and 
a\{l + 1) = pi- These and (I128p imply 


^1 Pdi (^a\) _ P‘1 pl—1 

C{l,Pd,oil) W-A 

The last two displays imply (jl38p for a(l) = 1. If / > a(l), let j > 1 be the position of I in 

o, i.e., I = a{j). In this case, the definition ()132p implies that the running products in the 

definitions of c* and c*^ have the same number of ratios and they are all equal except for 
the terms, which is {pi — X)/{pi — Pa{j-i)) for the former and {pi — X)/{pi — Pi-i) for the 
latter. ai has one more element than a, therefore, the ratio of the parity terms is again — 1 ; 
these imply 

Cg _ Pi - Pi-i 
Cgi Pi - Pa{j-1) ' 

On the other hand, Z e a, j > 1, oi = a u — 1} and (113311 imply a*{l) = p*{a{j — 1)), 

a*{I + 1) = pi, al{l) = pi-i, and a*{I + 1) = pi and therefore 

C'(^) Pdl Q^gj ) _ Pi Pl—1 

C{l,Pd,a*) Pl-Pa{j-1)' 

The last two displays once again imply (11381) when a{j) = I with j > 1- 
For a cz M, the definition (|133l) implies 

= «au{d}(^ + 1); 

we have already shown e H, then. Lemma [ 6 ] and the last display imply e Hi 

for I ^ a. Then by (I13ip e Hi for each loop on the vertex a u {d} of Ga, i.e., the last 

part of Definition 15.21 is also satisfied. This finishes the proof of the proposition. □ 

Proposition 6.4. 

d = 1,2, 3,..., d, are dB-determined Y-harmonic functions. 

Proof. That Zi^ is Y -harmonic follows from Propositions 16.31 and 15.21 The components of 
a cz { 1 , 2,3,..., d — 1} and /3J = pd are all between 0 and 1. This and Proposition 15.31 
imply that Zi^ are all 5i?-determined. □ 

With definition (11391) we can rewrite (jl35l) as 

-( n 2 (140) 

^ d 

for y e dB. 

Proposition 6.5. 

Py{T < oo) = 2 f n 1 ( 441 ) 

d=l \l=d+l - yd J 

for y e B. 
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Display (I163D below in subsection 18.11 shows the right side of (I14ip for d = 2. 

Proof. Let 1 e denote the vector with all components equal to 1. The decomposition of 
Gd into the single vertex {d} and d < d implies that the right side of (I14ip equals 




d-l 


d-l 


2 2 <m,al),y\ + Y, 

d=l aGV^d d=l 



- A 

fJ-i - ^J-d 


Kiy) 


for y e dB ; (11401) implies 

= [{PdA),y] 

which, for y e dB, equals 1. Thus, we see that the right side of (I14ip equals 1 on dB. 
Proposition 16.41 savs that the same function is 5i?-determined and is T-harmonic. Then its 
restriction to B must be indeed equal to y ^ Pyi'^ < oo); 2 / ^ -B, which is the unique function 
with those properties. □ 

7 Convergence - initial condition set for X 

In Section [3] we have proved a convergence result which takes as input the initial position of 
the Y process. One can also provide, as is done in the LD analysis, an initial position to the 
X"' process as X"'(0) = \nx\ for a hxed x e with XiieAG ^ ^ prove a convergence 
result in this setting. The initial position X"'(0) = \nx\ implies that probabilities such as 
those in (13 ip will all decay to 0 and therefore convergence to 0 no longer suffices to argue that 
a probability is negligible, we will now compare LD decay rates of the probabilities which 
appear in the convergence analysis. In the current literature only some of these rates have 
been computed in any generality. We believe that, at least for the exit boundary dAn, all of 
the necessary rates can be computed but forms a nontrivial task and will require an article 
of its own. Thus, instead of treating the general case, for the purposes of this paper, we will 
conhne ourselves to the case of two tandem queues in our convergence analysis when the 
initial position is given as X"'(0) = \nx\. 

In the rest of the section X will refer to the two dimensional tandem walk. The possible 
increments of X are vq = (0,1), = (—1,1) and V2 = (0, —1) with probabilities p{vq) = A, 

p{vi) = yi and p{v2) = ^2- For this model the stability condition (fTTp becomes A < yi,y2- 
On {x : x{l) = 0} [{x : x{2) = 0}] the increment (—1,1) [(0,-1)] is replaced with (0,0). 
For the present proof it will be more convenient to cast the limit in terms of the original 
coordinates of the X process. The Y process in the coordinate space of the X process is 
X = Tn{Y). X is the same process as X except that it is constrained only at the boundary 
^ 2 . 


Xk+I = Xfc + Tr{Xk,Ik) 

^k+i = d<-k+ 7ri(Xfc, /fc). 

We will assume that X and X start from the same initial position 


Xo = Xo 

and whenever we specify an initial position below it will be for both processes. 
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As before, r = inf{A: : Yk e dB} = e Tn = {k : Xi{k) + X2{k) = dAn} 

(by definition, X hits dAn exactly when Y hits dB); the subscript of P will denote initial 
position, i.e, Px{t < oo) equals P{t < oo) when Xq = Xq = x. 

Proposition 7.1. Let X and X be as above and assume A < ^i,/r 2 - For x e 0 < 
x(l) + x(2) < l,x(l) > 0 set Xn = \nx\. Then 

\Pxn{'^n < Tq) — PxnjT < 00)| fl42l 

{rn < To) 


decays exponentially in n. 

The proof will require several supporting results on cti = inf {A: : X^ e di} and 


(Ti ,2 = inf {A: : A: > ui, 6 ^ 2 }, 

di ,2 = inf{A; : k > fJi,Xfc(l) = -Xk{2)}. 


Proposition 7.2. 


for k ^ (Ti, 2 - 


Xfc(l)+Xfc(2) =Xfe(l)+Xfc(2) 


(143) 


Proof. 


Xk = Xk 


(144) 


for k ^ ai implies (11431) for A; ^ ui. If ui = < 71^2 then we are done. Otherwise X„.^{2) = 
Xo-j (2) > 0 and Xk{2) > 0 for cji < A; < let ai = ui < 1^2 < ■ ■ ■ < vk < cri^2 be the times 
when X hits di before hitting d2. The definitions of X and X imply that these are the only 
times when the increments of X and X differ: X^.+i — X^,. = 0 and Xj,^.+i — X{nj) = (—1,1) 
if = (—1,1); otherwise both differences equal luj. This and (|144p imply 


Xk — Xk = (,k ■ (— 1 ; 1 ) 


(145) 


for k ^ iTi ^2 where 


K 






and • denotes scalar multiplication. Summing the components of both sides of (|145p gives 
(fTi3l) . □ 


Define 

Tn = {cTi < <71,2 <Tn< Tq}. 

r„ is one particular way for {r^ < tq} to occur. In the next proposition we find an upper- 
bound on its probability in terms of 

7 = -(log(pi) V log(p2)) 

Proposition 7.3. For any e > 0 there is N > 0 such that if n > N 

PxA^n)^e-^(^-^\ (146) 


where Xn = [nxj and x e M+, a;(l) -I- x{2) < 1. 
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The proof will use the following definitions. 


Ha{q) = - log I 2 + S 1 ’ 

.is{0,l,2}—a {iGa} 


(147) 


where (•, •) denotes the inner product in For x e set 


b(x) = {i : x{i) = 0}. 


We will write H rather than 

Let us show the gradient operator on smooth functions on with V. The works MM 
use a smooth subsolution of 


^b(x)(VF(x)) =0 


(148) 


to find an lowerbound on the decay rate of the second moment of IS estimators for the 
probability Px^i^n < "To)- V is said to be a subsolution of (I148p if F7j,(x)(V^(3^)) ^ 0- The 
event T^ consists of three stages: the process first hits di, then 82 and then finally hits dAn 
without hitting 0. To handle this, we will use a function (s,a;) ^ 14(s,a;), with two variables; 
for the X variable we will substitute the scaled position of the X process, and the discrete 
variable s e {0,1,2} is for keeping track of which of the above three stages the process is 
in; V will be a subsolution in the x variable and continuous in (s,a;) (when (s,x) is thought 
of as a point on the manifold A4 consisting of three copies of (one for each stage); the 
zeroth glued to the first along di and the first to the second along ^ 2 ) and therefore one can 
think of V as three subsolutions (one for each stage) glued together along the boundaries 
of the state space of X where transitions between the stages occur. We will call a function 
(s,a;) ^ V{s,x) with the above properties a subsolution of (11481) on the manifold M.. 

Define 


y/(x) =<r„x> + 27 -( 3 -z)e, V^^^ = /\Vf, (149) 


i=0 


where 

ro = (0,0),ri = -7(l,0),r2 = -7(1,1)- 

The subsolution for stage j will be a smoothed version of As in mu, we will need 

to vary £ with n in the convergence argument; for this reason, e will appear as the third 
parameter of the constructed subsolution. The details are as follows. 

The subsolution for the zeroth stage is 1/(0, x,e) = 7 — 3e, Vl/(0, •) = 0 and it 

trivially satisfies (I148p and is therefore a subsolution. 

Define the smoothing kernel 


mix) 


■p^vix/S), r]{x) = l{|x|s:i}(|a;p 


1),M 



To construct the subsolution for the first and the second stages we will mollify 
j = 1, 2, with 7]: 

Vij,x,£)=( VX^{y)r]c^^{x-y)dy, (150) 

Jr2 

and C 2 is chosen so that 


1/(1,a;,e) = 1/(2,a:,e) 


(151) 
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for X e 82 and 


V{l,x,e) = V{0,x,e) (152) 

for X e di (this is possible since V{j, e, x) ^ as C 2 ^ 0 and all of the involved functions 
are affine; see m page 38] on how to compute C 2 explicitly). That V{j,-,e), j = 1,2 are 
subsolutions follow the concavity of Ha and the choices of the gradients rf, for details we 
refer the reader to m Lemma 2.3.2]; a direct computation gives 


^ £3 

dxidxj e ’ 


(153) 


j = 1,2, for a constant C 3 > 0 (again, the proof of |14l Lemma 2.3.2] gives the details of this 
computation). 

The construction above implies 


V{2,x,e) < 0,x e {x : a:(l) + x{2) = 1}. 
Now on to the proof of Proposition 17.31 
Proof. P(0, ■,£) maps to a constant and thus 

(VW{x),Vi} = W{x + Vi) - W{x) 


(154) 


(155) 


if IT = T(0, •,£). For IT = T(j, •,£), j = 1,2, Taylor’s formula and (I153p give 


VW(x),-v, 

n 





< 


£3 

ne 


(156) 


We will allow e to depend on n so that ^ 0 and n£n ^ 00 . Define Sk = l{o-i>fe} + Ifo-i 2 >fe}! 
Mo = 1 and 


Mk+i = MkBxp \-n[V [ S'fc+i, 


-V[Sk,—.£r 

n 


- 1 


{n>(Ti} 


C 3 

n£„ 


That V{j,-,£n), j = 0,1,2 are subsolutions of (I148p . the relations (I155p . (11561) (11521) and 
()151l) imply that M is a supermartingale (11561) and 1 (11551) allow us to replace gradients in 
(|148p and (11471) with finite differences and (|151l) and (I152p preserve the supermartingale 
property of M as S' passes from 0 to 1 and from 1 to 2). This and M ^ 0 imply (see [6l 
Theorem 7.6]) 




exp ^Sfe+i, 


V(Sk,—,e, 

n 


^{n>tTi} 


C 3 

nSn 


^ 1 , 


where ro,n = a tq. Restrict the expectation on the left to lr„ and replace l{n>o-i} with 1 
to make the expectation smaller: 


E.^ 


T0,ti 


Lr„e 


-T 0 ,ti 


exp 


J^vls. 


k=l 


X, 


’fe+ 1 : 


fe+1 



^ 1 . 


Over r„, X first hits di and then 62 and finally dA^. Furthermore, the sum inside the 
expectation is telescoping across this whole trajectory; these imply that the last inequality 
reduces to 


IE:r 


lr„e '*-""'°''‘exp(-n(T(2,X^ £„) -T(0,Xo,£n)) 


^ 1 . 
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To,n = Tn on and therefore on the same set ^to,„ ^ This, ^(0, •,£„) = 7 — 3e,i, (I154p 
and the previous inequality give 

^g-n(7-3£„)^ (157) 




lr„e 


~'^0,n 


Now suppose that the statement of Theorem 17.31 is not true, i.e., there exists e > 0 and a 
sequence Uk such that 

(158) 

for all k. Let us pass to this subsequence and drop the subscript k. m Theorem A. 1.1] 
implies that one can choose C 4 > 0 so that P{To^n > nci) ^ for n large. Then 


E„ 


lr„e 


--^3_T-n 

TO n 


^ [lr„e 

^ [lr„l{ro.„=£nC4}] 


P(£'i n E2) ^ T’(Ei) — P{E^ for any two events E\ and E2] this and the previous line 
imply 


— C 3 C 4 

^ g Tien {Pxni^n) Pxn{'^0,n > ^C4)) 
^ ^e-’T’(7-£) _ g-(7+l)>T.^ ^ 


By assumption nSn 00 which implies C3C4^jn£n 0; this and the last inequality say 


E, 


lr„e 


-X0,r, 


cannot decay at an exponential rate faster than 7 —e, but this contradicts 
(|157l) because £n —>■ 0. Then, there cannot be e > 0 and a sequence {n^} for which (I158|) 
holds and this implies the statement of Proposition 17.31 □ 


Define = log(p2)(l,l) and V{x) = (-log(/ 9 i) +<ri,a:» a (-log(p2) +<r 3 ,x», for 
X e 


Proposition 7.4. 

lim-logPa;„(rn < To) = V{x) 

n^co n 

for X e 0 < x(l) + x{2) < 1 and Xn = \nx\. 

The omitted proof is a one step version of the argument used in the proof of Proposition 
O and uses a mollification of V as the subsolution. 

Proposition 7.5. For any e > 0 there is N > Q such that if n > N 

PxifTl < (Ti ^2 < T < 00 ) ^ g-”'(7-e) (159) 

where Xn = [nxj and x e x(l) + x(2) < 1 . 

Proof. Write 

Px[(Ji < cri,2 < r < oo) = Px{cFi < cJi,2 < d-1,2 < r < oo) + Px[(Ti < (Ti, 2 < r < 0-1,2). 

The definitions of X and X imply tq ^ 0 - 1 , 2 . Then, if a sample path uj satisfies o-i(a;) < 
o'i, 2 (w) < r(a;) < 01 , 2 , it must also satisfy o-i(a;) < 0 - 1 , 2 (w) < Tn{<w) < To(a;). This and 
Proposition 17.31 imply that there is an N such that 

T’x„(o-i < 0 - 1,2 < r < d- 1 , 2 ) ^ 
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for n> N. On the other hand, Proposition 16.51 and the Markov property of X imply 

< 0 - 1,2 < ^ 1,2 < r < oo) ^ 

for some constant C 5 > 0. These imply (I159p . □ 

Proof of Proposition \n\ Decompose Px{Tn < tq) and Px{f < co) as follows: 

Pxr,{Tn < To) = Pxr,{rn < fTi < Tq) + (o'! < Tn ^ (Tp2 A Tq) (160) 

+ Pxr^if^l < 0-1,2 <rn< To) 

Pxr,{r < OO) = Px^{t < O-i) + Pa:„(o-l < T < fJp2) + PxA^^l < 0-1,2 < T < OO). (161) 

By definition X and X are identical nntil they hit di] therefore {r^ < ai] = {t < ai] and 

PxA'^n < o-i) = Px^T < ui). (162) 

The processes X and X begin to differ after they hit di; but Proposition 17.21 savs that 
the snms of their components remain equal before time 0 - 1 ^ 2 ; this implies r = on ^ < 71^2 
and therefore 

-Pa;„(o-l < T ^ 0 - 1 , 2 ) = PxS^^l < Tn ^ 0-1,2 A Tq) 

This (|162p and the decompositions (|160ll and (|161ll imply 

\Pxr,{jn < To) - Px^(t < 00)| = lPx^(cTl < 0-1,2 < Tn < Tq) - Px^io^l < 0-1,2 < T < 00)| 

By Propositions 17.31 and 17.51 for e > 0 arbitrarily small the right side of the last equality is 
bounded above by when n is large. On the other hand, Proposition 17.41 says for 

eo > 0 arbitrarily small Px„(Tn < tq) ^ e-^('n+<^o) fQj. large where 71 = V(x) < 7 . Choose 
e and cq to satisfy 

7 - 7i > e + eo. 

These imply that for cg = (e + eo) + 71 — 7 < 0 

\Pxn{Tn < Tq) — Px^jT < 00 )| ^ 

\Pxn{Tn < T))! 

when n is large; this is what we have set out to prove. □ 

It is possible to generalize Proposition 17.II in many directions. In particular, one expects 
it to hold for any tandem walk of finite dimension with the same exit boundary; the proof 
will almost be identical but requires a generalization of Proposition 17.41 which, we believe, 
will involve the same ideas given in its proof. We leave this task to a future work. 

There is a clear correspondence between the structures which appear in the LD analysis 
(and the subsolution approach to IS estimation) of pn and those involved in the methods 
developed in this paper. This connection is best expressed in the following equation (in the 
context of two tandem walk just studied): For q = {qi, 52 ) e set /3 = and a = 
then 

H[q) = -log(p(/3,a)), 

where p is the characteristic polynomial defined in ([96]). A similar relation exists between 
H2 and p 2 . In the LD analysis H and Hi appear as two of the Hamiltonians of the limit 
deterministic continuous time control problem; the gradient of the limit value function lies on 
their zero level sets. In our approach, the counterpart of H is the characteristic polynomial 
p; its 1 -level set H is the starting point of our definition of harmonic systems whose solutions 
give harmonic functions for the limit unstable constrained random walk Y of our analysis. 
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7.1 How to combine multiple approximations 

We have seen with Proposition 17.11 that ^Vn (r < oo) approximates Px^ijn < tq), Xn = [nxj 
extremely well (i.e., exponentially decaying relative error) for all x e ^4 = {x e M^,x(l) > 
0,x(l) + x{2) < 1} when n is large. In general this will not be true and to get a good 
approximation across all A we will have to use the transformation as well as T^. Thus, 
for general X, we will have to construct two limit processes and will be as above 

and Y^ will be the limit of = T^{X). In words, we obtain from X, by moving the 
origin to ( 0 , n) via an affine change of coordinates and removing the constraint on the second 
coordinate. In d, dimensions we will have d possible limit processes, one for each corner of 
dAn- A key question is how to decide for which of these limit processes Py^ix < oo) best 
approximates Px^i^n < t'o)- For the exit boundary we think that taking the maximum 
of the alternatives will suffice. We believe that the proof of this claim will involve arguments 
similar to those given above. We hope to provide its details in a future work, starting with 
the two dimensional case treated in this section. An example is given in subsection 18.21 


8 Examples 

We look at three examples: two tandem walk, general two dimensional walk and d-tandem 
walk. We have used Octave [7] for the numerical computations in this section and the rest 
of the paper. 


8.1 Two dimensional tandem walk 


Let us begin with the two tandem walk for which (11411) becomes 

PJr < 00) . pf (163) 

This gives the following approximation for Pxijn < tq): 

h -2 — A ^n-(3;(l)+3;(2))^3:(2) fJ-2 — ^ ^n-{x{l)+x{2)) ^x{2) (164) 

Proposition 17.11 says that for x e and Xn = [nxj, the relative error 


If(Xn) - PxA'Tn < To)| 

Pxr, {Xn < To) 

decays exponentially in n. Let us see numerically how well this approximation works. Set 
fii = 0.4, = 0.5, A = 0.1 and n = 60. In two dimensions, one can quickly compute 

Pxn{xn < To) by numerically iterating (p5]l and using the boundary conditions VsAn = 1 
and 1L(0) = 0; we will call the result of this computation “exact.” Because both / and 
Pxijn < To) decay exponentially in n, it is visually simpler to compare 

Vn = --log Px{xn < To), and Wn = --log/(x). (165) 

n n 

The first graph in Figure [10] are the level curves of Wn of Vn] they all completely overlap 
except for the hrst one along the x(2) axis. The second graph shows the relative error 
{Wn — Vn)lVn] we See that it appears to be zero except for a narrow layer around 0 where it 
is bounded by 0 . 02 . 
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0.02 



Figure 10: On the left: level curves of Vn (thin blue) and Wn (thick red); on the right: the 
graph of {Wn - Vn)lWn 

For X = (1,0), the exact value for the probability Px{tqo < tq) is 1.1285 • 10“^^ and 
the approximate value given by (I164p equals 1.2037 • 10“^^. Slightly away from the origin 
these quantities quickly converge to each other. For example, Px{Teo < tq) = 4.8364 • 10“^^, 
f{x) = 4.8148 • 10-35 for x = (2,0) and Px{m < tq) = 7.8888 • lO'^i, f{x) = 7.8885 • lO’^i 
for X = (9, 0). 

For X e and g : Z^ ^ M let T>g denote the discrete gradient of g: 




P(l) P(l) 

Figure 11: T>Vn{x) and PWn{x) , x = (5, •); on the right the same functions for x = (•, 1) 

{Vg){x) = ( 5 r(x + (l,0)) - g{x),g{x + {Q,l) - g{x)). 

The large deviations analysis of Vn suggests that nVVn approximately equals (pi,0) in a 
region around the x(l) axis and {p2,P2) elsewhere. These discrete gradients play a key role 
in the importance sampling estimation of the probability Px{Tn < tq). Since [H], it has been 
of interest to the author to understand how nVVn{x) transitions from (/9i,0) to {p2,P2) as x 
moves from the x(l)-axis to the interior of An- The approximation oi Vn hy f also explains 
how this transition takes place. As an example let us graph the values of these gradients 
over the line x(l) = 5 (any x(l) value slightly away from 0 will give similar results). The left 
panel of Figure fTTl shows the discrete gradients Wn and VWn along this line; they overlap. 
The right panel of the same figure shows the same gradients over the line x{2) = 1. 
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8.2 General two dimensional walk 


Now let us consider the two dimensional network with the transition matrix 

/ 0 0.15 0.1\ 

p = i 0.2 0 0.1 . (166) 

VO.24 0.06 0 / 

For this example, we will need to use both T^, i = 1,2, to get good approximations of 
Px{Tn < tq) across all of An. These two transformations will give us two limit unstable 
processes and T^. The first will give good approximations along 82 and the second along 
di. To combine their results into a single function, we will take their maximum. 

The limit processes and Y^ have the following dynamics: 

= + * = 1,2, 

where J* = J of (|2ip with i = 1,2. 

Remark 6. equals after we exchange the node labels (i.e., node 1 becomes 2 and 
2 becomes 1) This allows one to use the same computer code to compute either of the 
approximations by reordering the elements of the matrix p. 


We want to compute 

Py{T^ < cc),Py{T^ < 00 ), 

where r® = inf {A: : Y^ e dB}. We no longer have explicit finite formulas for these as we 
did in the tandem case. We will instead use a linear combination (a superposition) of basis 
functions of subsection sza to approximate the function mapping dB to 1; the same linear 
combination of the Balayage of the basis functions (for which we have explicit formulas) will 
provide an approximation for the probabilities we seek. One way to do this is as follows (we 
give the details for Py{T^ < 00 ), the procedure is identical for Pyir'^ < 00 ) ). As a first order 
approximation we use 


An = 


C(r,a(r, 1)) 


hr = [(r, 1),-] 


C{r,l) 

C{r,OL{r,l)) 


[(r,Q:(r,l)),-]. 


By Proposition 15.11 Aq is a harmonic function of Y^. 

For the p of (I166p . /3i(l) = r = 0.42373 and Q:(r, 1) = 0.48123. Then, by Proposition 
14.101 Aq is dB determined. These imply 


Py{T^ < 00 ) — AqI ^ Py{T^ < 00 ) max I Ao(y) — 1 

yedB 


(167) 


Set 

c, ^ 

C(r,Q:(r, 1)) 

Aq — 1 = C 7 [(r, Q:(r, 1)), •] is geometrically decreasing on dB and therefore it takes its greatest 
value at n = 0 where it equals, for the present example, 3.8418 < 4. This and (I167p imply 


^Ao < Py{T^ < 00 ) < Aq. 

0 

Thus, even with a single T-harmonic pair of log-linear functions, we are able to approximate 
Py{T^ < 00 ) up to a constant term. To improve, approximate 


C7[r,a{r,l))r] 
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by a superposition of harmonic pairs of subsection 14.71 as follows. Consider the vector 
b = C 7 ([(r, Q:(r, 1)), (y, y)],y e {0,1, 2,..., ii'}) e If one thinks of the restriction of 

C 7 [(r, a{r, 1)), •] to dB as a sequence, one truncates it to its first K+1 components to get &; for 
K large (for the present example we take K = 11), the remaining tail of C 7 [(r, a(r, 1)), -Jigs 
(its components from the K + 2”''^ on) will be almost 0. What we want to do is to construct 
basis vectors Vi,i = 0,1,2, for by truncating in the same way the restrictions 

to dB of ii' + 1 log-linear y-harmonic functions and write 6 as a linear combination of 
the members of this basis. To construct our first basis element take the harmonic function 
[(/3(ri),ri), •] of Proposition 112] and define vq = ([(/?(ri), n), (y, y)], y e {0,1, 2,..., K}). 
This gives us a vector in to complete it to a basis for we need K more vectors. 

Set aj = i?e "+i, where Re (ri, 1) is to be specified and consider the T-harmonic functions 
^j3i(aj) of Proposition 15. 11 We would like all of these to be di?-determined, for which 

\Pi{aj)\,\aj\,\a{/3i{aj),aj)\ < 1 (168) 

suffice; the second of these is satisfied by definition. The sufficient conditions we have 
derived for the rest, listed as Proposition 14.131 don’t cover the parameter values of the 
present example (([56]) and p(0,2) = 0 fail). Then, what we will do is to compute them 
explicitly (using ([8^ for a{l3i{aj),aj) and (f53]l for and verify directly that (I168p 

holds. Figure [12] shows the results of these calculations for R = 0.7 and K = 11 and indeed 
we see that |/3i(aj)|, \ oL[f5i{aj), aj)\ < 1 holds for all j. Thus, by Proposition l4.10l all 




5 10 15 20 25 30 35 40 


Figure 12: /3i(Q;j)’s (shown with x’s) and Q:(/3i(aj), aj)’s (shown with ’o’s) on the C-plane; 

the graph of the error A defined in 


are di?-determined. Dehne 


K+l 


Vj = (V(«i)(y’y)’y = o>i>2,...,ii') e c 

Dehne the change of basis matrix B to consist of rows vq, directly evaluating its 

determinant shows that B is invertable (this determinant is a polynomial in aj and /3j, this 
can be used to show that, perhaps after perturbing aj, we can always assume B invertable). 
Dehne the coefficient vector 

= bB \ 


By dehnition. 


K 


^1 = V’(0)[(/3(ri),ri),-] 2 V’(j)V(a^) 

1=1 
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equals c^\{r, cx{r, 1)), •] over the set {(y, y), y = 0,1, 2, K} cz dB. That |Q:(r, 1)| <1, (I168|) 
and 0 < ri < 1 imply that 

A(y) = I Ai(y, y) - C 7 [(r, cx{r, 1)), (y, y)]| ^ 0 (169) 

exponentially as y ^ oo. Then one can explicitly find a bounded interval \^K + 2, it''] in 
which A(y), y s Z+, takes its maximum value. For K = 11 and for the parameter values of 
the current example, this difference takes its maximum value at y = 12 (see the right panel 
of Figure [12]) where it equals 0.00796 < 0.008. These imply 

0.992(Ao + Ai) < Py{T^ < oo) < 1.008(Ao + Ai), y e B. 

Set gi = Ao +Ai. Using exactly the same ideas we construct a function 52 ( 2 /) approximating 
Pyij'^ < 00 ). gi and 52 give two possible approximate values for Px{Tn < tq): gi{T^{x)) and 
g 2 {T^{x))] as pointed out in subsection 17.11 one expects 

f{x) = max(yi(r^(x)),y 2 (F^(a:))) 

to be the best approximation for Pxijn < tq) that one can construct using gi and 52 - As 
in the previous section, we compare / and Px{Tn < tq) in the logarithmic scale. Define 
Wn{x) = -^logf{x) {Vn{x) is, as before, Vn{x) = -^\ogPx{Tn < To)). Figure m shows the 
graph of iVn — Wn)lVn for the present case: qualitatively it looks similar to the right panel 



Figure 13: Relative error for the nontandem two dimensional example 

of Figure [101 the relative error is near zero across A„, except for a short boundary layer 
along the a;(2)-axis; where it is bounded by 0.03. One difference is the slight perturbation 
from 0 of the relative error on dA^ which comes from the approximation error depicted in 
the right panel of Figure [T^ 

The level curves of Vn, —^\oggi{T^{x)) are shown on the left panel of Figure [TH those 
of Vn and — ^ log 52 (7]?(^)) given on the right panel. It is clear from these graphs that, 
indeed, as discussed above, gi{T^[x)) approximates Vn well away from d*. 

Finally, suppose we are given a function h on dA„. The algorithm above can also be used 
to approximate [/i(T^(Y^l))l{T-i<ooj], i = 1,2, and hence E^, [h{Xx^)l[r„<To}]- We leave 
the analysis of this approximation to future work. 
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Figure 14: Two approximations given by two limit processes 

8.3 Tandem walk in higher dimensions 

Take a four dimensional tandem system with rates, for example, 

A = 1/18,= 3/18,/r2 = 7/18,//s = 2/18,/i4 = 5/18. 

Let f{y) denote the right side of (|14ip . As before, define V/ = — log{Px{Tn < To))/n and 
Wn = — log f{Tn{x))/n. The level curves of Vn and Wn and the graph of the relative error 
{y _ W)IV for X = (0,0, z,j), i,j ^ n, are shown in Figure fTCl once again, qualitatively, 
these graphs show results similar to those observed for the earlier examples: almost zero 
relative error across the domain selected, except for a boundary layer along the x(4)-axis, 
where the relative error is bounded by 0.05. 




Figure 15: Level curves and relative error in four dimensions 

As our last example, consider the 14-tandem queues with parameter values shown in 
Figure [T6l 

For n = 60, An contains 60^^/14! = 8.99 X 10^^ states which makes impractical an exact 
calculation via iterating (1351) . On the other hand, (11411) has 2^'^ — 1 = 16383 summands and 
can be quickly calculated. Define Wn as before. Its graph over {x : x(4) -|-x(14) = Q0,x{j) = 
0 ,j A 4,14} is depicted in Figure [T71 For a finer approximation of T/i q,--- ,o)('^n < tq) 
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Figure 16: The service rates (blue) and the arrival rate (red) for a 14-dimensional tandem 
Jackson network 




Figure 17: The graph of Wn over {x : x(4) -I- x(14) = 60, x(j) = 0,j ^ 4,14} 


we use importance sampling based on Wn- With 12000 samples, IS gives the estimate 
7.53 X 10“^*^ with an estimated 95% confidence interval [6.57,8.48] x 10“^° (rounded to two 
signihcant figures). The value given by our approximation (11411) for the same probability is 
/((1,0,-- - ,0)) = 1.77 X 10 which is approximately 1/4*^ of the estimate given by IS. The 
LD estimate of the same probability is (A/minl^^(/ij))®^ = 4.15 x 10“^^. The discrepancy 
between IS and (I14ip quickly disappears as a;(l) increases. For example, for a:(l) = 4, IS 
gives 2.47 x 10”^® and (11411) gives 2.32 x 10“^®. 

9 Conclusion 

The foregoing analysis points to a number of future directions for research. We state some 
of them here. 
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9.1 Constrained diffusions with drift and elliptic equations with Neumann 
boundary conditions 

Diffusion processes are weak limits of random walks. Thus, the results of the previous 
sections can be used to compute/approximate Balayage and exit probabilities of constrained 
unstable diffusions. We give an example demonstrating this possibility. 

For a, 6 > 0 let X be the the constrained diffusion on M x M_|_ with infinitesimal generator 
L defined as 

/ - Lf, Lf = <V/, ((2a + 6), (a - 6))> + j) ’ 

where denotes the Hessian operator, mapping / to its matrix of second order partial 
derivatives. On {x : x(2) = 0} X is pushed up to remain in M x M_|_ (the precise dehnition 
involves the Skorokhod map, see, e.g., m)- o, 6 > 0 implies that, starting from B = {x : 
x(l) > x(2)}, X has positive probability of never hitting dB = {x : x(l) = x(2)}. Let r be 
the hrst time X hits {x : a;(l) = x{2)}. Proposition 16.51 for d = 2 suggests 

P^(t < OO) = g-{‘^o,+b)3x{2) 

* a — b 

_ aP^^-3(2a+b)x(l) ^ 3 , e H. (170) 

a — b 

One can check directly that the right side of the last display satishes 

LV = 0, {VP, (0,1)) = 0,x e 62 - 

This and a verification argument similar to the proof of Proposition 14.101 will imply (|170l) . 
Almost the same argument for general d gives an explicit formula for the solution of the 
d dimensional version of the above elliptic equation on M x with 2^~^ — 1 Neumann 
boundary conditions on d{M x 

9.2 Solutions to perturbed nonlinear PDE 

As indicated in the introduction, classical large deviations analysis leads (at least for Markov 
processes) to a deterministic first order HJB equation. To improve the approximation pro¬ 
vided by the solution of this first order PDE, one can add nonlinear second order perturbation 
terms to it m- We expect the ideas of the paper to bear on the task of computing approx¬ 
imate solutions of the perturbed second order nonlinear PDE related to the probabilities 
treated in the present paper. 

9.3 Extension to other processes and domains 

In the foregoing sections, we have computed approximations to the Balayage operator and 
exit probabilities of a class of constrained random walks in two stages: 1) use an affine 
change of coordinates to move the origin to a point on the exit boundary and take limits; 
as a result, some of the constraints in the prelimit process disappear and one obtains as a 
limit process an unstable constrained random walk; 2) find a class of basis functions on the 
exit boundary on which the Balayage operator of the limit process has a simple action; then 
try to approximate any other function on the exit boundary with linear combinations of the 
functions in the basis class. The type of problem we have studied here is of the following 
form: there is a process X with a certain law of large number limit which takes X away 
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from a boundary dAn towards a stable point or a region; tq is the first time the process gets 
into this stable region. We are interested in the probability Pijn < tq) and the associated 
Balayage operator. We expect the first step to be applicable to a range of problems that fit 
into this scenario. The second stage obviously depends on the particular dynamics of the 
original process and the geometry of the exit boundary. It remains to be explored for which 
processes and boundaries it is possible to construct classes of simple basis functions. Even 
very simple changes in the dynamics or in the boundary geometry from those covered in this 
paper may lead to different types of basis functions. We hope to treat the tandem walk case 
with a separate boundary in an upcoming work. 

9.4 Harmonic functions with polynomial terms 

Let us confine ourselves, for the purposes of this brief comment, to tandem queues. If 
= /U 2 , (IhTll no longer holds and indeed (II63p is not well defined (because of the /ii — ^2 
in the denominator of the first ratio). One way to remedy this is to replace the right side of 
(116311 by the limit of the same expression as /ii —> /i 2 ) which gives 

Py{T < (») = + ^(y(l) - ^ (171) 

fl 

where /ri = /r 2 = A* and p = X/p. Similarly, in three dimensions one gets, for example, for 
Pi = P2 = P3 

Py{T < 00 ) = ((^ + y(3)co^ + co^ y{l) + 1^ , 

where cq = {p — X)/p and y{l) = p(l) — (y(2) + p(3)). Similar limits can be computed 
explicitly for the cases pi = P 2 Ps, Pi = P 3 p 2 and pi ¥= P 2 = P 3 - Generalization of 
these results to higher dimensions and more general topologies remain for future work. 

9.5 The boundary layers of Vn and subsolution based importance sampling 
algorithms 

The works mu develop IS algorithms based on subsolutions of (I148p to estimate pn- 
These works and others which followed them express most of their functions in a law of 
large numbers scale (as we do in Sections [3] and [7]) . We will express everything in unsealed 
coordinates in the discussion below. Again, to be brief, we will limit ourselves to formal 
comments on two tandem queues. Details and generalizations remain for future work. 

For certain values of system parameters, Vn of (jI65p may manifest a boundary layer, 
where its discrete gradient sees a rapid change near the boundaries of its state space. This 
happens when pi = p 2 for the case of two tandem queues with a boundary layer along the 
x(l) axis. Here is one interpretation of the subsolution approach of [SJ [H] in the present 
context (i.e., two tandem queues and pi = p 2 )'- the discrete gradient of the large deviation 
value function 

Vn{x) = -Uog 

approximates the discrete gradient of Vn well everywhere except along a boundary layer along 
the x(l) axis; the subsolutions constructed in [5l[T3] are perturbations of Vn which attempt 
to approximate the discrete gradient of Vn also in this boundary layer. The subsolutions 
constructed in these works involve a parameter that satisfy 

TiCn * QC, (^n * 0 (f72) 
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and have boundary layers of constant width parallel to the x(l) axis; e„ determines the 
width of the boundary layer. Although O Theorem 3.8] says that (|172l) suffices for the IS 
algorithm defined by the subsolutions to be asymptotically optimal, m subsection 2.3.3] 
observes that a good performance of the algorithm for finite n (accurate estimation results 
with bounded estimated relative error) requires a finer specihcation of and hence of the 
size of the approximating boundary layer. How to measure the size of the boundary layer 
of Vn (and use this information to specify e^) more precisely has remained an open problem 
since the inception of the subsolution approach. With (I17ip we now see that the correct way 
to add a boundary layer to Vn (so that it has a boundary layer mimicking that of Vn) is to 
perturb it to 

Wn{x) = --log (p^-VW+V^)) + - (x(l) + . (173) 

n V ) 

Proposition 17.II and the numerical example of subsection 18.11 suggest that the boundary layer 
of Wn matches that of Vn as n increases. 

The definitions (|149p and (|15np give an alternative construction of smooth subsolutions 
with explicit boundary layers (the region where = ri). In contrast, the boundary layer 

of Wn is expressed implicitly in its definition (117311 . Let us now try to quantify explicitly the 
size and the shape of the boundary layer of Wn and thus that of Vn- 
Define 

W{y) = - log + ^(2/(1) - y(2))p"(')) . 

Then Wn{x) = ^W{Tn{x)). It suffices to calculate the boundary layer of W] this we can 
transform by Tn to get that of Wn- We will specify the boundary layer by its boundary 
I : R+ ^ M+; the layer will be defined as {y : y{2) ^ Z(y(l))}. The defining property of 
the layer is that it is the region where the gradient of W rapidly changes. Away from the 
boundary {y(2) = 0} W behaves like the linear function {y{l) — y{2)) log(/3) whose directional 
derivative V(i^i)lT in the direction (1,1) is zero. And indeed the same is true for W itself 
on {y : y{l) = y{2)}, i.e., 

dW dW 

V(M)(1P) ^ ^(2/) + ^(2/) = 0 ,ye{y: y{l) = y{2)}- 

Furthermore, for hxed 2/(1), ^( 1 , 1 )^[u) is decreasing in y{2) and the above display implies 
that this directional derivative hits zero on di- We will define l{y{l)) as the point where the 
value of V(ip)W is half of its value at (y(l),0), i.e., 

Z(2/(l)) = 2/(2)* (174) 

where 2/(2)* is the unique solution of 

fe(l) - !,(2)*) (1 + 1^9(1)) - . 

I is increasing 2 /( 1 ): this implies that when transformed by Tn it defines a boundary layer for 
Wn (and hence for Vn) that narrows down as it extends toward the point (n, 0). In contrast, 
the subsolutions developed in [5l [H] have boundary layers of constant size, i.e., parallel to 
the a;(l) axis. The graph of x(l) ^ Z((n — x(l))), x(l) e [0, n], and the level sets of Vn for 
n = 40, A = 0.2, and y = 0.4 are shown in Figure [T8l 
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Figure 18: The contours of Vn for n = 40, A = 0.2 , ^ = 0.4 and its boundary layer 
computed using I of (I174p 
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